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PREFACE 

This  document  is  designed  as  a companion  document  to  AFWL-TR-75-165, 
SIGANAL:  A MODULAR  SIGNAL  ANALYSIS  PROGRAM  developed  by  Mr  Ramon  A.  Tenorio, 


Computational  Services  Division,  Air  Force  Weapons  Laboratory.  Mr  Tenorio 's 
programs  greatly  streamlined  and  consequently  simplified  hardware  data  analysis 
procedures  at  the  Laboratory.  This  software  immeasurably  increased  the  capa- 
bil«ty  for  application  of  sophisticated  analysis  procedures  to  our  data  analysis 
Aided  by  this  new  capability  and  strongly  encouraged  by  Mr  Tenorio,  the  author 
set  himself  to  the  task  of  demonstrating  the  utility  and  interpretation  of 
modern  random  data  analysis  capability  afforded  by  the  new  software. 

My  thanks  to  Tony  Tenorio  for  his  strong  encouragement  and  his  uniquely 
organized,  efficiently  and  easily  used  analysis  package.  His  efforts  in  pro- 
gramming several  of  the  specialized  routines  used  in  the  simulation  analysis  is 
gratefully  acknowledged. 

Many  hours  of  interesting  and  enligntening  discussions  with  Dr  Paul  Merritt 
have  clarified  my  understanding  of  signal  coherence  analysis.  I thank  him  for 
his  insight. 
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SECTION  I 

INTRODUCTION  AND  GUIDE 


This  guide  Is  especially  tailored  for  use  as  c companion  document  to  SIGANAL: 
A Modular  Signal  Analysis  Program  (reference  1).  SIGANAL  Is  a user's  guide  to  a 
comprehensive  package  of  software  routines  for  signal  analysis  applications. 

These  routines,  coded  In  FORTRAN  IV,  process  a user  specified  sequence  of  data 
manipulation  routines  which  operate  on  digitized  data.  Presently  coded  routines 
facilitate  calculation  of  signal  statistics  Including  means,  standard  deviations, 
signal  value  probability  density  and  distribution  functions,  Fourier  transforms 
(discrete),  power  and  cross-spectral  densities,  coherence  functions  and  linear 
system  transfer  functions.  This  guide  demonstrates  the  usefulness  to  which  these 
routines  can  be  applied  In  extracting  pertinent  signal  information. 

1 PURPOSE  OF  SIGNAL  ANALYSIS 

Signals  embody  and  represent  characteristics  of  the  systems  which  generate 
them.  Signals  generated  by  sensors,  for  example,  convey  basic  information  about 
the  process  being  measured.  The  signal  information,  conveyed  typically  as  an  elec- 
tric voltage  or  current,  represents  some  physical  quantity  such  as  position,  veloc- 
ity, acceleration,  force,  pressure,  temperature,  luminosity,  brightness,  field 
strength,  fluid  flow,  volume,  weight,  length,  density,  charge,  or  any  of  innumer- 
able other  quantities.  In  some  cases,  the  analyst  is  concerned  only  with  static 
measurements  wherein  a single  number  adequately  represents  the  desired  quantity. 
That  the  gravitational  acceleration  on  the  earth's  surface  Is  9.81  meters/sec2  is 
sufficiently  accurate  for  almost  everyone  except  the  experimental  and  the  gravity 
wave  physicists.  We  shall  not  be  concerned  with  these  constant  signals.  We  are 
Interested  In  fluctuating  signals  representing  time-varying  and  dynamic  system 
variables.  Fourier  analysis  and  other  transform  techniques  have  been  developed 
to  represent  these  signals  In  a frequency  domain  basis  particularly  suitable  for 
calculating  the  response  of  linear  dynamics  systems  or  analyzing  requirements  of 
communication  channels.  Techniques  of  stochastic  signal  analysis  reviewed  in  this 
. guide  are  suitable  for  these  purposes;  however,  their  true  utility  resides  in  in- 

terpretation of  randomly  fluctuating  signals  which  arise  in  phenomena  where  future 
, events  are  not  deterministically  predictable  from  past  events.  These  signals 

occur  In  statistical  mechanics,  quantum  mechanics,  turbulence  and  noise  theories. 
This  guide  defines  parametrlzatlons  of  stochastic  process  and  descriptive  param- 
eters which  may  be  measured  or  estimated. 
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We  adopt  the  philosophy  that  power-spectral  densities  (PSD's)  are  the  most 
convenient  format  for  characterizing  stochastic  processes.  We  emphasize  this 
approach  throughout.  It  has  evolved  from  our  experience  that  linear  dynamic  sys- 
tems described  by  differential  equations  have  Input-output  properties  that  are 
most  readily  characterized  and  viewed  as  frequency  domain  transfer  functions 
rather  than  convolution  Integrals.  Frequency  domain  characterizations  of  both 
system  forcing  functions  and  system  transfer  functions  provide  ready  visualiza- 
tion of  output  signal  characteristics  (In  the  frequency  domain).  By  our  emphasis 
on  the  power-spectral  density,  we  devote  little  attention  to  discussion  of  auto- 
correlation functions  other  than  to  define  them  and  to  derive  their  relationship 
with  power  spectra.  Our  preference  Is  not  universal. 

2 USER'S  GUIDE  TO  SIGNAL  ANALYSIS 


Concepts  essential  to  the  understanding  and  Interpretation  of  practical 
stochastic  signal  analysis  procedures  and  results  are  presented  In  Sections  II 
and  III.  Section  IV  concludes  these  developments  with  the  presentation.  Interpre- 
tation and  comment  on  the  use  of  SIGANAL  routines  to  analyze  random  data  synthe- 
sized by  a numerically  simulated  dynamical  system  excited  by  random  disturbance 
phenomena . 

Basic  mathematical  foundations  defining  and  relating  properties  of  stochastic 
processes  are  presented  in  Section  II.  Introductory  probability  theory  is  briefly 
developed  and  expanded  to  define  stochastic  processes.  Probability  concepts  are 
generalized  to  stochastic  process  characterizations  by  autocorrelation  functions 
and  power  spectral  densities.  Attention  Is  restricted  to  ergodlc  processes  which 
most  practical  stochastic  processes  approximate.  The  mathematical  basis  for  calcu- 
lating PSD's,  cross-spectral  densities,  coherence  functions,  and  transfer  function 
estimates  is  developed. 

In  Section  III,  the  practical  aspects  of  calculating  stochastic  process  char- 
acterizations are  reviewed.  These  aspects  Include  data  sampling,  finite  duration 
observation  Intervals,  and  confidence  bounds.  The  Fast  Fourier  Transform  technique 
for  efficiently  calculating  signal  spectra  Is  Introduced  and  specific  calculation 
formulas  presented.  These  topics  review  technique  as  implemented  in  the  SIGANAL 
code,  and  accent  the  proper  application  of  analysis  procedures  and  interpretation. 

Section  IV  Is  devoted  to  the  development  of  a simple  numerical  signal  analysis 
problem  which  exemplifies  the  basic  Interpretation  principles  and  augments  our 
Intuition  as  to  what  to  expect  in  the  way  of  results.  A first  order  control 
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system  excited  by  error  sensor  noise  and  an  exponentially  correlated  disturbance 
process  Is  simulated.  Theoretical  power  spectra  and  coherence  functions  are 
calculated  from  the  known  system  transfer  functions  and  the  white  noise  distur- 
bance processes.  These  "expected"  results  are  compared  to  the  quantities  calcu- 
lated with  the  signal  analysis  routines. 

Readers  having  no  prior  familiarity  with  data  analysis  should  proceed  sequen- 
tially from  Section  II  through  Section  IV  and  should  augment  his  reading  with 
material  In  the  references.  Those  more  knowledgeable  with  signal  analysis  pro- 
cedures should  scan  Section  II  to  familiarize  themselves  with  notation  (see  also 
the  nomenclature)  and  then  may  study  the  numerical  example,  referring  back  to 
Section  III  as  required  for  reference. 
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SECTION  II 

PROBABILITY  AND  STOCHASTIC  PROCESSES 


1 INTRODUCTION 

A detailed  exposition  of  probability,  stochastic  processes,  and  random  data 
analysis  as  presented  In  the  references  (2,  3,  4)  Is  beyond  the  scope  and  intent 
of  this  guide.  We  do  present  the  basic  definitions  and  concepts  required  for  an 
introductory  understanding  of  these  subjects,  particularly  In  their  application 
toward  analyzing  random  effects  observed  In  measurements.  First,  the  notions  of 
random  effects  are  introduced  to  explain  phenomena  which  deterministic  system 
concepts  cannot  predict.  So  motivated,  we  characterize  stochastic  processes  and 
how  they  Interact  with  our  systems. 

2 BASIC  CONCEPTS:  RANDOMNESS 

The  basic  notion  of  random  data  or  random  signals  stems  from  our  Inability 
to  deterministically  characterize  signals  or  phenomena  we  observe.  To  determin- 
istically characterize  a signal,  we  mean  that  we  can  specify  expllclty  signal 
values  as  a function  of  time  In  units  of  voltage,  amperage,  pressure  or  some 
other  appropriate  unit.  For  example,  we  know  that  a voltage  V of  10  volts  im- 
pressed upon  a 1 pf  capacitor  In  parallel  with  a 1 Mn  resistor  will  decay  as 
V(t)  = 10  e-t  once  the  voltage  source  Is  removed.  The  voltage  has  been  expressed 
as  a well-known  mathematical  function  of  time.  Another  Important  aspect  of  the 
example  Is  that  we  obtain  the  same  result  whenever  the  experiment  is  repeated. 

Not  all  experiments  are  precisely  repeatable.  The  radio  Interference  from 
an  automotive  Ignition  system  does  not  duplicate  with  precisely  the  same  output 
radio  signal.  Other  experiments  Involving  fluid  turbulence,  electronic  noise  and 
molecular  or  atomic  collisions,  for  example,  result  in  phenomena  which  cannot  be 
precisely  predicted.  These  phenomena  are  random  because  they  cannot  be  precisely 
specified.  In  subsequent  sections,  statistical  characterizations  of  these  pro- 
cesses shall  be  defined.  These  characterizations  shall  serve  as  a basis  for  de- 
scribing the  response  of  systems  having  such  random  excitations. 

Development  of  analysis  definitions  and  techniques  Is  required  for  an  under- 
standing of  the  response  of  dynamic  systems  to  random  excitations.  Application 
of  these  techniques  Is  often  essential  for  proper  system  design  and  system  per- 
formance analysis.  Just  as  a system  designer  can  compute  precisely  the  control 
system  response  to  a deterministic  Input  (e.g.,  step  function,  sinusoid. 
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exponential)  we  wish  to  be  able  to  compute  the  response,  on  the  average,  to  ran- 
dom Inputs.  Further,  we  must  know  the  average  variation  by  which  one  particular 
response  may  differ  from  the  average  response.  These  data  allow  specification  of 
system  sensors  and  actuators  so  that  they  will  not  be  saturated  by  the  noises  and 
disturbances  affecting  the  system.  Statistical  characterizations  of  random  dis- 
turbances and  their  effect  on  dynamical  systems  also  allow  the  designer  to  com- 
pute system  accuracy  or  performance  bounds  for  those  cases  In  which  random  dis- 
turbances are  the  fundamental  system  limitations.  With  this  motivation,  we  begin 
with  a discussion  of  elementary  probability  theory.  A comprehensive  treatment  Is 
available  In  Papoulls  (reference  4). 

An  experiment  whose  outcome  cannot  be  predicted  Is  said  to  be  random.  The 
result  of  a coin  toss  Is  therefore  random  since  we  cannot  predict  In  advance 
whether  the  result  will  be  HEADS  or  TAILS.  There  are  some  important  aspects  of 
this  experiment,  however,  that  can  be  described.  For  a fair  coin,  we  would  ex- 
pect that  the  result  HEADS  Is  just  as  likely  as  the  result  TAILS.  Intuitively 
then,  we  would  expect  that  If  the  coin  were  tossed  N times,  HEADS  would  come  up 
approximately  N/2  times  (l.e.,  NH  « N/2)  and  TAILS  would  come  up  N-j. « N/2  times. 

The  ratio  N^/N  Is  called  the  relative  frequency  of  the  outcome  HEADS  in  the  coin 
toss  experiment.  We  define  the  probabl 1 1 ty  of  the  event  HEADS  as  its  relative 
frequency  In  the  limit  as  the  number  of  coin  tosses  approaches  infinity.  Notation- 
ally, 

P ({HEADS})  = 11m  Nu/N 

| H 

N-*» 

(2-1) 

Certain  properties  of  the  coin  toss  experiment  are  common  to  all  random  ex- 
periments. These  properties  are  the  basis  of  probability  theory.  Each  experiment 
has  a set  of  mutually  exclusive  results  or  outcomes.  To  each  outcome  we  assign  a 
positive  number  (possibly  zero)  representing  the  probability  (relative  frequency) 
that  the  outcome  will  occur  If  the  experiment  Is  repeated.  Sets  of  outcome  we 
call  events  and  to  each  event  we  assign  a probability  equal  to  the  sum  of  the 
component  outcome  probabilities.  The  event  containing  all  possible  outcomes  has 
probability  P - 1,  These  concepts  are  fully  developed  in  reference  4,  Chapter  2, 
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3 RANDOM  VARIABLES 

To  each  outcome  of  an  experiment  we  may  also  assign  an  arbitrary  number 
called  a random  number.  The  functional  relationship  between  the  number  x and  the 
events  ? Is  denoted  x(e)  and  the  function  x is  called  a random  variable,  provided 
that  It  satisfies  certain  general  conditions.  Basically  these  conditions  entail 
the  requirement  that  the  set  {x^y}  corresponds  to  a set  of  experimental  out- 
comes (l.e.,  an  event).  In  our  coin  toss  experiment,  for  example,  we  arbitrarily 
define  a random  variable  as  the  rule  which  assigns  value  1 to  the  outcome  HEADS 
and  value  0 to  the  outcome  TAILS.  Me  could  have  just  as  easily  assigned  the  value 
1/2  to  HEADS  and  1/2  to  TAILS,  as  the  particular  values  assigned  are  not  of  Impor- 
tance. A key  concept  of  the  random  variable  Is  that  there  exists  a "one-to-one" 
correspondence  between  experimental  events  and  sets  of  random  numbers.  Thus 
{x  * 1}  corresponds  to  the  events  HEADS  and  {x  = 0},  or  more  generally  {x  < 1}, 
corresponds  to  the  event  TAILS.  In  each  case  the  random  number  associated  with 
each  outcome  is  contained  In  the  given  set  of  x values. 

Since  each  set  of  random  numbers  Is  associated  with  an  event  of  the  under- 
lying probability  experiment,  we  associate  with  each  set  the  properties  of  the 
underlying  events.  Most  Importantly,  we  associate  with  each  set  of  random  num- 
bers the  probability  P of  the  associated  event.  The  sets  {x^y}  play  a particu- 
larly important  role.  Any  set  of  random  variables  associated  with  an  event  can 
be  expressed  in  terms  of  the  basic  sets  {x^y}  related  by  set  operators  union. 
Intersection  and  complement.  Properties  of  general  sets  are  then  obtained  in 
terms  of  properties  of  the  basic  component  sets.  Important  properties  of  these 
random  variable  sets  are  presented. 

The  probability  distribution  function  is  defined  as  the  probability  of  the 
event  {x-^y}  and  written 


F(y)  = P({x^y}) 


(2-2) 


The  probability  distribution  function  for  the  random  variable  defined  on  the  coin 
tossing  experiment  Is  plotted  In  figure  1. 

The  probability  density  function,  pdf.  Is  defined  as  the  derivative  of  the 
probability  distribution  function.  The  density  function  for  the  coin  toss  experi- 
ment Is  plotted  In  figure  2.  Notatlonally, 

f 0)  - 3y  f(y) 
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Figure  1.  Coin  Toss  Experiment  Probability 
Distribution  Function. 


Figure  2.  Coin  Toss  Experiment  Probability 
Density  Function. 


The  random  variable  associated  with  the  coin  tossing  experiment  Is  a discrete 
random  variable  since  all  events  correspond  to  at  most  a countable  number  of  spe- 
cific random  variable  values.  Discrete  random  variables  have  pdf's  which  are 
collections  of  Impulse  functions. 

More  generally  an  experiment  may  have  a continuum  of  possible  outcomes.  The 
probability  density  function  for  these  random  variables  Is  a continuous  function 
except  at  possibly  a countable  number  of  points.  The  pdf  for  a uniformly  dis- 
tributed random  variable  Is  Illustrated  In  figure  3.  The  corresponding  probabil- 
ity distribution  function  Is  plotted  In  figure  4.  This  random  variable  Is  uni- 
formly distributed  since  each  value  Is  equally  likely  to  occur.  Suppose  we  wish 
to  know  the  probability  that  {yj  < x^y2}.  By  definition  we  know  that 
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Figure  3.  Uniformly  Distributed  Random 


Variable  - pdf. 


Figure  4.  Uniformly  Distributed  Random  Variable  - 
Probability  Distribution  Function. 

F(yi)  = P Ux^yO) 

F(y2)  * P ({x^y2}) 

Because  the  event  {x^y2}  Includes  the  event  {x^y^  and  since  {yx 
equals  {x  < y2}  less  {x  < yj},  it  follows  that 

P({yi x ^ y2})  = P({x^y2})  - P({x^yil) 


F(y2)  - F(yi) 
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By  Integration  of  equation  (2-3)  we  have  that 


■/. 


F(y)  * / f(w)dw 


(2-7) 


Integral  properties  and  equations  (2-6)  and  (2-7)  allow  us  to  express  equation 
(2-6)  as 


P({yi 


. f 

Jyi 


f(w)dw 


(2-8) 


and  for  the  case  that  y2  * yx  +Ay,  Ay  being  a differential  quantity,  we  have 


P({yi^x^yx  +Ay})  = f(yx)  Ay 


(2-9) 


Thus  f(y)  represents  a differential  probability  since  It  can  also  be  defined  as 
the  limit 


f(yx)  =>  11m 

Ay  -*•  0 


P({yi  < x < yx  + Ay}) 
Ay 


(2-10) 


The  properties  of  the  probability  distribution  and  density  functions  ex- 
pressed In  equations  (2-4)  through  (2-10)  are  true  in  general.  Additional  proper- 

/ 

ties  are  presented  In  reference  3,  Chapter  2,  and  reference  4,  Chapter  4.  For  the  uniform 
density,  f(y)  Is  a constant  for  0 < y < 1 so  that  by  equation  (2-10)  the  probabil- 
ity of  any  differential  Interval  of  random  numbers  Is  the  same.  That  Is,  each 
differential  Interval  Is  equally  likely  or  has  probability  zero  If  the  Interval 
falls  outside  [0,1] . 

The  normal,  chi-square,  binominal,  beta,  F-dlstributlon  and  student's  t- 
dlstrlbutlon  are  among  those  encountered  In  practice.  We  limit  our  attention  to 
the  normal  and  the  chi-square  distributions.  These  distributions  play  an  Impor- 
tant role  In  random  data  analysis.  The  normal  distribution  Is  particularly  Impor- 
tant In  random  variable  theory  since  any  random  variable  which  Is  the  sum  of  K 
Identically  distributed  Independent  random  variables  has  a density  function  that 
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approaches  a normal  distribution  as  K Increases.  This  result  Is  guaranteed  by 
the  central  limit  theorem  (reference  5).  The  distribution  and  density  functions 
for  a normally  distributed  random  variable  are  plotted  In  figures  5 and  6,  re- 
spectively. The  normal  density  function  Is  written: 


f(y) 


1 


; 


exp  [-1/2  (*^)2] 


(2-11) 


Figure  5.  Normal  Random  Variable  Probability 
Distribution  Function. 


; 


Figure  6.  Normal  Random  Variable  Probability 
Density  Function. 


The  density  function  Is  completely  specified  by  the  parameters  u and  o.  The 
plausibility  of  the  central  limit  theorem  Is  demonstrated  in  figure  7 in  which 
density  functions  for  the  function  WK  equal  to  the  normalized  sum  of  uniformly 
distributed  Independent  random  variables  are  plotted. 
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Figure  7.  Demonstration  of  the  Central  Limit  Theorem. 


VK  1=1 


(2-12) 


Also  plotted  In  figure  7 Is  a normal  density  function.  Observe  that  the  density 
function  for  approaches  the  normal  density  as  K increases. 

4 MEANS  AND  MOMENTS 

Except  for  a few  special  cases.  It  Is  usually  Inconvenient  or  mathematically 
cumbersome  to  characterize  random  variables  by  specifying  their  probability  dens- 
ity or  probability  distribution  functions.  These  functions  can  be  difficult  to 
specify  whenever  they  are  not  expressible  In  terms  of  known  or  tabulated  functions 
as  Is  done  In  equation  (2-11).  The  pdf  totally  characterizes  the  random  variable. 
An  alternative  approach  to  completely  characterizing  the  random  variable  through 
Its  probability  density  function  Is  to  examine  exactly  what  properties  of  the  ran- 
dom variable  we  are  Interested  In.  Then  we  partially  characterize  the  random 
variable  In  terms  of  parameterlzatlons  of  these  useful  properties.  This  approach 
Is  exemplified  through  the  following  definitions  and  Illustrative  examples. 
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Two  quantities  are  typically  of  most  Interest  In  the  investigation  of  random 
phenomena.  A primary  quantity  of  Interest  is  the  mean  or  average  value  of  the 
random  variable.  That  Is,  if  we  were  to  observe  Independent  repetitions  of  the 
same  experiment,  the  numerical  average  of  these  observations  is  the  single  most 
Important  quantity  for  us  to  compute.  The  second  most  important  quantity  is  the 
dispersion  or  spread  of  the  observations  about  the  mean  value.  These  concepts 
are  quantified  as  follows. 

We  define  the  expectation  operator,  denoted  £[•],  as 


<£H 


[*]  f(y)  dy 


(2-13) 


where  the  term  within  the  brackets  Is  a specified  function  of  the  random  variable 
y.  We  Interpret  the  expectation  operator  as  evaluating  the  average  value  of  the 
function  upon  which  It  operates.  The  average,  of  course,  is  with  respect  to  the 
random  variable  and  the  underlying  probability  space  upon  which  It  Is  defined. 

Let  us  examine  this  Important  concept  from  another  point  of  view.  The  average 
value  of  any  quantity  Is  simply  the  sum  (Integral)  over  all  possible  values  that 
the  quantity  may  assume  of  a product  formed  as  a value  multiplied  by  the  proba- 
bility that  the  particular  value  Is  assumed.  Notationally 


Average  (g) 


P({g  <x^  g + dg>) 


(2-14) 


Application  of  equation  (2-10)  allows  us  to  re-express  equation  (2-14)  as 


Average  (g) 


9 f(g)  dg 


(2-15) 


We  now  realize  that  equation  (2-13)  Is  simply  a variation  of  equation  (2-15).  The 
equations  differ  In  the  following  respect.  The  bracketed  quantity  In  equation 
(2-13)  Is  a function,  say  g(y),  of  y.  Particular  values  of  g ma^  be  obtained  by 
greatly  different  values  of  y.  Rather  than  combining  the  probabilities  of  any  y 
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which  will  yield  the  specified  value  g to  obtain  an  equation  of  the  form  equation 
(2-15),  we  recognize  that  we  can  evaluate  the  probability  density  for  each  y sep- 
arately and  simply  Integrate  over  all  possible  y values  to  obtain  equation  (2-13). 
We  Illustrate  this  with  an  example. 

EXAMPLE  2-1 


Suppose  we  wish  to  calculate  the  average  of  the  absolute  value  of  the  random 
variable  y which  Is  uniformly  distributed  In  the  Interval  [-1/2, 1/2] . We  have 
that: 


f(y) 


1 , -1/2  — y — 1/2 
0 , otherwise 


g(y)  ■ |y| 


By  direct  application  of  equation  (2-13), 

€[\y\]  * / |y|  f(y)  dy 

*1/2 


■L 

■f 


|y|  dy  = 1/4 


1/2 

A pdf  for  g Is  obtained  by  first  calculating  the  probability  of  the  random  vari- 
able g. 


P({x  ^ g^  x + dx})  * 


P({x  ^ y ^ x+  dx})  + P ( {x  -y  x+  dx»  , x^O 
0,  otherwise 


Then  by  application  of  equation  (2-10)  we  compute  that: 


f(g)  * 2 f (y)  , g^O 

0 , otherwise 


Finally,  evaluation  of  equation  (2-15)  gives 
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AVERAGE  (g) 


2gdg  = 1/4 


We  have  obtained  the  same  result  by  application  of  equation  (2-13)  or  equation 
(2-15).  Often  In  practice  we  shall  find  that  equation  (2-13)  is  easier  to  evalu- 
ate, generally  because  f(g)  may  not  be  readily  expressed  or  as  readily  calculated 
as  It  Is  In  this  example.** 

The  expectation  operator  provides  a basis  for  defining  two  sets  of  quanti- 
ties. The  nth  moment  of  the  random  variable  y is  defined  as 


[y"] 


f(y)  dy 


The  first  moment  is  the  mean  or  average  value  of  the  process  and  is  denoted  by 


Hy.  The  nth  central  moment  is  defined  as 


^[(y'wy)n] s f (y-*y)n  f(y>dy 

✓ —00 


The  first  central  moment  is  always  equal  to  zero.  It  is  of  no  interest  to  us. 
The  second  central  moment  is  called  the  variance  of  the  distribution  and  is  de- 
noted by  Oy2.  The  square  root  of  the  variance  is  called  the  standard  deviation 
of  the  process.  Both  the  variance  and  the  standard  deviation  characterize  the 
dispersion  of  the  distribution  about  its  mean  value.  The  larger  the  variance  or 
standard  deviation,  the  more  probable  It  Is  that  we  will  observe  values  of  the 
random  variable  outside  a fixed  Interval  about  the  mean. 

Higher  order  moments  usually  are  not  of  as  much  importance  to  us  as  are  the 
mean  and  standard  deviation  of  a process.  Generally,  however,  the  odd  central 
moments  (l.e.,  n equal  to  an  odd  Integer)  characterize  the  skewness  or  asymmetry 
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of  the  random  variable  probability  density  function.  Density  functions  which 
satisfy  the  symmetry  condition 


f(uy+a)  = f (yy-a)  , for  any  a 


have  Identically  zero  odd  central  moments.  The  even  central  moments,  like  the 
second  central  moment,  characterize  the  dispersion  of  the  distribution.  We 
Illustrate  the  calculation  of  mean  and  variance  with  the  following  example. 

EXAMPLE  2-2 

Calculate  the  mean  and  variance  of  a normally  distributed  random  variable 
with  the  density  expressed  In  equation  (2-11).  Explicitly 


l.  y ^2?  6XP  ^"V2  ^ 


Observe  that  the  density  function  Is  symmetric  about  y so  that 


(y-u)  f(y)  dy  = 0 


and  since  the  Integral  of  the  density  function  is  unity 


u / f (y)  dy  = y 


y f (y)dy  = yy 


Thus  the  quantity  y specified  In  the  normal  density  function  Is  the  mean  value  of 
the  distribution.  To  obtain  the  variance  we  note  that: 


f(y)dy  * 1 


and  multiplying  both  sides  by  a we  have  expllclty 
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(y-u)2  f(y)  dy  = o2 

Thus  the  quantity  a2  specified  In  the  normal  density  function  is  the  variance  of 
the  distribution.  The  normal  density  function  Is  completely  specified  In  terms 
of  the  mean  u and  the  variance  a2. 

Plots  of  normal  density  functions  having  the  same  mean  and  different  variances 
are  presented  In  figure  8.  Observe  that  a larger  variance  of  the  density  func- 
tion corresponds  to  a larger  dispersion  or  spread  of  the  density  function  about 
the  mean./-* 

For  a given  probability  experiment,  there  is  no  unique  random  variable  having 
real  values  assigned  to  each  outcome  or  event  of  the  experiment.  Infinitely  many 
random  variables  may  be  defined  upon  the  same  experiment.  Individually  these  vari- 
ables display  the  properties  discussed  in  section  II. 3 and  II. 4 of  this  report. 
Together,  two  or  more  random  variables  defined  upon  the  same  experiment  have  addi- 
tional properties  which  we  shall  find  useful  In  our  data  analysis.  Properties  of 
two  or  more  random  variables  are  defined  in  the  next  section. 

5 PROPERTIES  OF  SEVERAL  RANDOM  VARIABLES 

Analogous  to  equation  (2-2)  we  define  the  joint  probability  distribution 
function  of  two  random  variables  as  the  probability  of  the  event  {x^  y and  ws=z} 
and  is  written 

F(y,z)  * P({x^y  and  ws=£z})  (2-19) 
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fi*2  y 


Figure  8.  Variance  Effects  of  the  Normal  Density  Function. 

Similarly,  the  joint  probability  density  function  Is  defined  as  the  partial  de- 
rivative of  the  joint  distribution  function  with  respect  to  each  variable.  Nota- 
tionally. 


= wF(^z) 


The  random  variables  y and  z are  Independent  If  and  only  if  the  joint  density 
function  can  be  expressed  In  the  form* 


f(y.z)  ■ fy(y)fz(z) 


*The  subscripts  In  equation  (2-21)  distinguish  the  two  density  functions  fy(y) 
for  the  random  variable  y and  f2(z)  for  the  random  variable  z.  The  subscripts 
have  no  other  significance. 
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Independence  of  two  random  variables  means,  practically,  that  knowledge  of 
one  of  the  variables  does  not  convey  any  Information  regarding  the  other. 

Joint  moments  and  central  moments  are  defined  analgous  to  equations  (2-16) 
and  (2-17)  for  single  random  variables.  Explicitly,  joint  moments  are  defined  as 
the  expectation  of  the  product  of  Integer  powers  of  the  random  variables  and  de- 
noted 


Joint  central  moments  are  defined: 


(2-22) 


Cki  "f[(y-»y)k(I-ui)*I 


(2-23) 


The  order  n of  these  moments  Is  defined  as  the  sum  of  the  subscripts.  For  example 
m12  Is  a third  order  moment;  Cu  Is  a second  order  central  moment.  The  second 
order  central  moment  Cn  Is  of  particular  Importance  In  data  analysis  and  Is 
called  the  covariance.  For  fixed  variances  of  the  component  random  variables, 
an  Increase  In  Cu  magnitude  corresponds  to  a greater  and  greater  linear  depen- 
dence between  the  two  random  variables.  The  linear  dependence  of  two  random  vari- 
ables Is  characterized  directly  by  the  correlation  coefficient  pyz  which  Is  simply 
Cu  normalized  by  the  product  of  the  standard  deviations  of  the  component  processes 


p 


yz 


*ii 


a a 
y z 


(2-24) 


♦The  definition  of  expectation  presented  In  equation  (2-13)  must  be  generalized 
to  the  case  of  two  random  variables  for  application  to  equation  (2-22)  above. 
Explicitly  the  expanded  definition  is 


H f(y,z)  dy  dz 


Expectation  means  the  average  with  respect  to  all  the  component  random  variables. 
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Two  random  variables  are  linearly  dependent  If  one  of  the  random  variables  can 
be  expressed  as  a linear  function  of  the  other  plus  a third  independent  random 
variable.  In  the  case  that  no  such  linear  function  exists,  the  two  random  vari- 
ables are  linearly  independent  or  uncorrelated.  Linearly  Independent  random 
variables  have  Identically  zero  correlation  coefficients.  Thus  from  equation 
(2-24)  we  conclude  that: 


c“  ’ ^[(y"“y)  (2'“y)] 


= £[yz]  - u u2  = 0 


(2-25) 


Often  equation  (2-25)  Is  used  to  define  uncorrelated  random  variables  as  those 
random  variables  satisfying  the  condition  that 


£[yz]  * uyuz 


(2-26) 


Independent  random  variables  are  always  uncorrelated  since  equation  (2-26)  fol- 
lows directly  from  equations  (2-21)  and  (2-22). 

We  explicitly  show  the  Important  relationships  between  correlation  coeffi- 
cients and  linear  dependence  of  random  variables  by  the  following  construction. 
Assume  that  random  variables  y and  z with  means  and  u2,  respectively,  and 
nonzero  variances  ay2  and  az2  and  a correlation  coefficient  pzy  are  given. 

Then  we  shall  show  that  there  exists  a random  variable  w such  that 


(2-27) 
(2-28) 

% * vy  " auz  (2-29) 

°w2  * °y2  t 1 " Pyz  ] (2-30) 


y-az,  1 .e.  y s az  + w 


■=  °JL 
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SK^z)  (w-"w)]  • 0 


(2-31) 


We  proceed  with  the  proof  of  these  properties  by  defining  a random  variable  x 
and  then  showing  that  x satisfies  all  the  properties  attributed  to  w In  equations 
(2-27)  through  (2-31)  above.  Because  functions  of  random  variables  are  also  ran- 
dom variables,  we  define  x as  the  difference  between  y and  a scalar  multiple  b of 


x ty  - bz 


(2-32) 


The  scalar  b Is  a constant  which  we  conveniently  choose  so  that  x satisfies  the 
desired  properties.  Clearly,  by  taking  expectations  of  equation  (2-32) 


u * u - bu 
x My 


(2-33) 


and  subtracting  this  result  from  equation  (2-32),  squaring  each  side,  and  again 
taking  expectations  one  obtains: 

a2  = <?{[(y-uy)  - b(z-uz)]2> 


- a-y  - 2bCn  + b2  o2 


(2-34) 


Now  we  explicitly  evaluate  the  covariance  of  x and  z and  equate  It  to  zero  by 
appropriate  choice  of  b. 

e[(z-»z)  K)]  ■#{(«-«,)  [(y-uy)  -!>(«-*,)]} 


Cu  -bo2  sit  0 


(2-35) 


from  which  we  obtain: 
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b = Cn/  o*  = ay/az  (Cxl/ayaz^ 


- 

= — 0 
“x  yz 


(2-36) 


Substitution  of  this  result  Into  equation  (2-34)  and  substitution  for  Cn  gives 


_ p2  a2 

x y Hyz  y 


(2-37) 


Thus  making  the  Identification  of  x with  w and  b with  a and  comparing  equations 
(2-32),  (2-36),  (2-33),  (2-37),  and  (2-35)  with  equations  (2-27)  through  (2-31), 
respectively,  we  have  completed  the  desired  proof.  These  observations  follow 
directly  from  equations  (2-27)  through  (2-31). 

(I)  as  the  correlation  between  y and  z Increases,  the  variance 
of  w,  a*  decreases.  For  p =1,  a*  = o 

W Jr  Z W 

(II)  |ptfJ  — 1 since  o*  must  be  nonnegative 

jr£  W 

(III)  < o*  and  equality  holds  only  If  p£_  = 0 

w y yz 

(1 v ) the  error  In  approximating  y by  a linear  function  of  z Is  pre- 

cisely the  variance  of  w;  the  correlation  coefficient  squared 
Is  precisely  the  power  In  y attributable  to  z normalized  by  the 
total  power  In  y.  Explicitly, 


n»2  (z-"Z)2J 


(2-38) 


Zero  correlation  between  random  variables  y and  z Implies  that 
the  linearity  coefficient  a of  equation  (2-27)  Is  Identically 
zero  and  that  w = y 

The  quantities  a,  p . and  are  uniquely  determined  by  the 
yz  w w 

means,  variances,  and  the  covariance  of  the  random  variables  y 
and  z 
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The  concepts  of  linear  dependence,  correlation  and  joint  random  variables 
and  their  application  to  a data  analysis  problem  are  demonstrated  In  the  next 
example. 

EXAMPLE  2-3 

Suppose  that  we  wish  to  estimate  the  value  of  an  unknown  resistance  R by 
simultaneously  measuring  the  current  I through  the  resistor  and  the  voltage  V 
across  It.  Then  using  Ohm's  law  we  obtain  the  resistance  as  R - V/I.  Further 
suppose  that  because  of  old  equipment  and  coarse  meter  dials  we  cannot  obtain  a 
sufficiently  accurate  value  for  R.  We  repeat  the  experiment  shown  plctorlally 
In  Figure  9 for  various  values  of  source  voltage  Vs  and  make  the  measurements 
tabulated  below.  From  these  values  plotted  In  figure  10  we  wish  to  extract  the 
"best  estimate"  of  the  unknown  resistance. 

V(volts)  .28  .29  .48  .69  .97  1.07  1.24  1.23  1.53  1.81 

Urna)  T22  7SA  “47  — 7TCD — TM — TTTS — 1725 — 1751 — VM 

We  formulate  a "least-squares"  solution  by  assuming  that  the  measured  voltage  V 
Is  the  sum  of  a linear  function  of  I and  a measurement  noise  term  n. 

We  write:  V * RI  + b 
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1 2 
CURRENT  (ma) 


Figure  10.  Mean  Square  Estimation  of  a 
Voltage  Current  Relationship. 

where  the  caret  above  V and  R denotes  that  they  are  estimated  quantities.  The 
constant  b denotes  instrument  bias.  Since  we  have  assumed  V differs  from  V by 
a noise  term,  we  have 

"1  " V1  - *i 

We  arbitrarily  choose  the  unknown  parameters  R and  b so  as  to  minimize  the  noise 
variance,  or  equivalently,  we  minimize  the  mean-square  error,  by  which  an  esti- 
mated quantity  V differs  from  the  actual  measurement  V, 

We  define  a performance  measure 


e*to  I h -K + b)]2 

1-1 


and  we  minimize  E by  optimum  choice  of  the  unknown  parameters  R and  b.  The  solu- 
tion of  the  unconstrained  minimization  problem  Is.most  readily  obtained  by  equat- 
ing the  partial  derivatives  of  E with  respect  to  R and  of  E with  respect  to  b 
both  equal  to  zero  and  solving  these  equations  simultaneously  for  the  unknowns  R 
and  b.  Executing  these  steps 
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3E/3R  = -2/10  £ (V1I1  - Rif  - blj  * 0 
1*1 


3E/3b  * -2/10  Y (V^l  - b)  a 0 
1*1 

and  solving  the  above  equations  simultaneously  yields  the  result 


1/10  Y.  Vr(  1/10  Y *1 )( 1/10  Y vi 

1*1  \ 1*1  / \ 1*1  > 


9.6Kfl 


-i/io  Y ji  1/10  E vi  M1/10  E 1/1°  E vi 


0.04  volts 


where 


10  / io  \2 

D = 1/10  Y !1  "(1/10  E M 
1*1  V 1*1  ] 


In  view  of  the  relative  frequency  Interpretation  of  probabilities  developed  In 
section  I I. 2,  we  may  re-lnterpret  the  above  equations  In  terms  of  the  expec- 
tations they  approximate.  We  have 
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i/io  £ 

1=1 


Vi*uv 


10 


1/10  £ - Of  + ,{ 


1=1 

10 


1/10 


1*1 


Thus  expressing  R and  b In  terms  of  these  approximations  gives 


" Cu/»f  * ^ PVI 


Vi  “ uicn 

b « -xl  -1—  * vv  - R », 

°f 


* 

c ; 


Had  we  known  the  statistics  of  V and  I we  could  have  simply  applied  equation 
(2-28)  to  have  calculated  R In  terms  of  ay,  aj,  and  pyj.  Typically,  however, 
these  quantities  are  not  known  a priori  but  must  be  approximated.  Equivalently, 
the  mean-square-error  procedure  developed  In  this  example  yields  the  best  answer 
given  the  available  data.** 

Thus  far  the  concepts  of  randomness  and  probability,  random  variables  and 
their  characterization  by  probability  density  or  distribution  functions  have  been 
developed.  The  mean  or  average  value  and  the  variance  or  dispersion  of  the  ran- 
dom variable  have  been  defined  and  shown  to  be  two  of  the  more  Important  charac- 
terizations. Covariance  and  the  correlation  coefficient  are  shown  to  be  Important 
properties  of  random  variables  whenever  functions  of  random  variables  are  of  In- 
terest. 
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6 STOCHASTIC  PROCESSES 

The  concept  of  random  variables  Introduced  In  section  II. 3 has  no  time  de- 
pendence associated  with  It.  This  limitation  Is  serious  when  one  Is  most  Inter- 
ested In  the  effect  of  random  phenomena  upon  systems  having  Input-output  rela- 
tionships described  by  differential  equations.  In  this  section,  the  concept  of 
stochastic  processes  Is  Introduced  to  allow  analysis  of  those  situations  In 
which  time  dependent  random  effects  excite  linear  dynamical  systems  or  corrupt 
a time  dependent  Information  signal  we  wish  to  examine. 

A time  function  x(t)  having  numerical  values  which  depend  not  only  upon 
the  value  of  Its  argument  t but  also  upon  the  outcome  of  probability  experiment 
Is  called  a stochastic  process.  For  example  the  function 


x(t) 


sin  wt  , If  HEAD'S 
cos  wt  , If  TAIL'S 


(2-39) 


where  HEADS  or  TAILS  Is  the  outcome  of  a coin  toss  experiment  as  In  section  II. 2 
Is  an  example  of  one  stochastic  process.  In  practical  cases  of  interest,  the 
underlying  probability  space  characterizes  the  size  distribution  of  atmospheric 
turbules  or  the  "shot"  rate  of  electronic  noise  sources  or  some  other  physical 
phenonema  having  outcomes  which  cannot  be  precisely  predicted  In  advance  but 
which  do  having  varying  effects  upon  our  particular  experiment.  Often,  the 
underlying  probability  space  Is  not  well  characterized,  known,  or  understood 
other  than  through  observation  of  the  resulting  signal  or  effect  on  our  experi- 
ment. The  observed  effect  Is  also  a stochastic  process  because  its  value  de- 
pends upon  a stochastic  process  excitation  which  in  turn  depends  upon  the  out- 
come of  the  underlying  probability  experiment  outcome.  The  observed  stochastic 
process  may  be,  for  example,  the  Intensity  fluctuations  I(t)  of  an  electromag- 
netic wave  propagated  through  a turbulent  atmosphere  or  the  voltage  variations 
of  an  electronic  system  excited  by  noise. 

Random  signal  analysis  techniques,  which  are  the  primary  topic  of  this 
guide,  are  In  essence  a collection  of  methods  for  characterizing  these  stochas- 
tic processes  and  the  transform  relationships  which  model  the  Interaction  of 
stochastic  signals  with  linear  dynamical  systems. 

Basically,  stochastic  processes  are  defined  on  probability  experiments  in 
the  same  manner  as  random  variables.  In  fact,  for  fixed  agrument  t stochastic 


34 


J 


f 

1 

' 


processes  degenerate  to  random  variables.  Characterizations  of  random  vari- 
ables are  adequate  to  model  certain  fixed  time  properties  of  stochastic  pro- 
cesses but  new  concepts,  definitions  and  characterizations  are  required  to 
parameterize  Important  additional  properties  of  the  stochastic  process. 

Recall  that  a stochastic  process  defines  a real  function  of  the  indepen- 
dent variable  t for  each  outcome  of  a probability  experiment.  The  collection 
of  all  possible  time  functions  Is  the  ensemble  of  the  stochastic  process.  Both 
the  member  of  the  ensemble  (or  equivalently  the  outcome  of  the  experiment)  and 
time  must  be  known  to  assign  a value  to  the  stochastic  process  function.  A 
stochastic  process  ensemble  Is  partially  depicted  in  figure  11  which  shows 


Figure  11.  A Stochastic  Process  Ensemble. 

explicitly  the  dual  dependence  of  the  stochastic  process.  The  properties  mean, 
variance,  nth  order  moments  and  central  moments,  joint  moments,  covariance,  and 
correlation  developed  In  sections  II. 4 and  II. 5 for  random  variables  are  equally 
applicable  to  a stochastic  process  at  each  fixed  time.  Consequently,  these 
properties  are  defined  as  functions  of  time  for  stochastic  processes.  Additional 
properties  of  the  stochastic  process  to  be  developed  in  this  section  are  general- 
izations of  the  basic  properties  that  have  already  been  introduced.  Stochastic 
processes  are  completely  characterized  by  joint  probability  distribution  func- 
tions of  the  form: 

F^Xj  , X2  ......  , X^  , t j , t2**  “ 

x(^l ) — xi » x(tz)  — *2».  • *x  ^ (2-40) 
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which  must  be  known  for  §21  finite  sets  of  {t^},  (reference  3,  p.  51).  Such 
functions  are  cumbersome  or  defy  specification  in  all  but  the  simplest  cases. 
Fortunately,  several  other  properties  of  stochastic  processes  — namely  sta- 
tlonarlty,  ergodlclty,  and  normality  — are  generally  valid  and  allow  a greatly 
simplified  characterization.  These  properties  are  precisely  defined  subsequent 
to  the  following  preliminary  definitions. 

7.  MEANS,  MOMENTS,  AND  ERGODICITY 

The  concepts  of  means  and  moments  developed  In  sections  I I. 4 and  I I. 5 for 
random  variables  are  explicitly  extended  In  this  section  to  stochastic  processes. 
First,  probability  density  functions  and  expectation  operators  are  defined  for 
stochastic  processes.  The  probability  distribution  function  equation  (2-40) 
for  stochastic  processes  Is  considered  for  the  case  consisting  of  precisely  a 
single  element  tj.  We  define  the  first  order  density.  In  analogy  with  equation 
(2-3),  as 


. 

I 


! 


\ 


(2-41) 

The  denstty  f(xlfti)  has  explicit  dependence  on  time  ti  by  virtue  of  the  time 
dependence  of  the  probability  distribution  function.  An  expectation  operator 
follows  from  equation  (2-41)  as  a straightforward  generalization  of  expectation 
for  a single  random  variable  as  defined  In  equation  (2-13). 


[•]  f(x,t)  dx 


(2-42) 


The  mean  Is  defined  as  the  expectation  of  the  process  Itself. 


Vx(t)  a£[x(t)]  =J  x(t)f (x,t)dx 


(2-43) 


The  mean  at  any  fixed  time  Is  simply  the  average  of  x over  the  elements  of  the 
ensemble. 
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Stationarlty  defines  the  property  that  the  probabilities  describing  our  pro- 
cess do  not  depend  upon  the  time  origin  or  the  time  frame  In  which  we  describe 
our  process.  For  example,  stationarlty  Implies  that  the  probability  density 
function  satisfies  the  equalities  Indicated  In  the  following  cases, 

ffXi.ti)  = f(x1,t1+x) 


f(Xl,X2,t1,t2)  i f(X!,X2,tl+T,t2+T) 


(2-48) 


which  must  be  valid  for  any  and  all  t.  Clearly,  then  f(x1,t1)  is  Independent  of 
tj  and  may  be  more  simply  written  f(xx).  Similarly,  f(xx,  x2,  tx,  t2)  Is  Inde- 
pendent of  the  specific  time  reference  used  and  therefore  depends  only  upon  the 
relative  time  difference  x = tx  - t2.  Thus  for  stationary  processes, 


f(xx,  x2,  tx,  t2)  = f(xx,  x2,  x) 


(2-49) 


Stationarlty  reduces  the  complexity  of  the  probabilistic  characterization. 
Stationary  processes  have  means  vx  which  are  time-invariant.  The  correlation 
and  covariance  functions,  like  the  stationary  second-order  density  functions, 
depend  only  upon  the  time  difference  x.  We  denote  these  functions  m (x)  and 

Ajr 

C (x)  for  the  stationary  process  correlation  and  covariance  functions,  respec- 
*«y 

tlvely.  We  shall  consider  only  stationary  processes  throughout  the  remainder  of 
the  guide. 

Expectations  have  been  defined  as  averages  over  the  underlying  probability 
experiment.  In  signal  analysis,  another  Important  average  Is  time-average. 
Time-average  mean,  correlation  and  covariance  quantities  are  defined  for  a par- 
ticular realization  x^  and  y^  of  the  stochastic  processes  as 


X1 


tlm 

T-*» 


x^(t)dt 


(2-50) 
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rV2 

R.v(t)  = Aim  1/T  / Xj(t+T)yj(t)dt 

y T-  J-T/2 


(2-51)* 


rT/z , 

Cxy(t)  = Aim  1/T J ^ [Xl  (t+x)  - xjfy^tj-yjdt  ^ 

provided  that  the  Indicated  limits  exist.  The  time-averages  defined  above  are 
random  variables  since  they  have  values  dependent  upon  the  particular  realiza- 
tion and  y^  of  the  underlying  stochastic  process  (reference  4,  p.  326). 
Assuming  process  statlonarlty  and  taking  expectations  of  the  above  quantities  we 
find  that 


Similarly, 


£[Rxy(r)]  = mxy(x) 


(2-53) 


(2-54) 


£[  Vt)]  = Cxy(x) 


(2-55) 


*In  the  literature  R (x)  Is  also  defined  as 

xy 


rT/2 

Rxy(x)  * Aim  1/T  / Xi(t)  y^t+x)  dt. 

T~  J-T/2 


The  definition  Is  simply  a matter  of  convention  and  does  not  substantively  change 
the  results.  The  reader  Is  cautioned,  however,  to  familiarize  himself  with  the 
particular  convention  adopted  by  another  author. 
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For  many  stochastic  processes,  the  random  variables  x,,  R (t)  and  C (t)  have 

• xy  xy 

zero  variance,  have  Identical  values  Independent  of  the  stochastic  process 
realization  and  used  In  their  time-average  calculation,  and  are  Identi- 
cally equal  to  the  corresponding  ensemble  average.  Should  these  three  proper- 
ties hold  for  all  average  values  of  a stochastic  process,  then  the  process  Is 
said  to  be  ergodlc.  An  Important  consequence  of  ergodiclty  Is  that  each  ele- 
ment or  the  stochastic  process  ensemble  Is  representative  of  the  ensemble  as 
a whole.  Ergodiclty  allows  us  to  compute  Important  characterizations  (mean, 
covariance,  etc.)  of  a stochastic  process  by  applying  time-average  calcula- 
tions to  a simple  observation  of  the  ensemble.  We  shall  consider  only  ergodlc 
stochastic  processes  throughout  the  remainder  of  this  guide.  We  state  without 
proof  that  most  of  the  stochastic  processes  encountered  In  engineering  practice 
are  ergodlc  or  nearly  ergodlc  In  the  sense  that  our  analysis  techniques  remain 
valid  and  useful. 

8 LINEAR  DYNAMIC  SYSTEMS 


A linear  dynamic  system  Is  a system  in  which  certain  dependent  (output) 
variables  satisfy  a given  linear  differential  equation  having  one  or  more  inde- 
pendent variables  or  forcing  functions.  We  shall  consider  the  general  case  In 
which  these  forcing  functions  are  stochastic  processes  and  derive  stochastic 
process  characterizations  of  the  system  output  In  terms  of  the  Input  and  the 
system  transfer  function  or  differential  equation.  Having  completed  this  char- 
acterization we  show  the  Importance  or  correlation  and  covariance  calculations 
In  estimating  dynamic  system  Input-output  relationships.  We  assume  the  reader 
has  a knowledge  of  differential  equations  and  Fourier  transforms  such  as  pre- 
sented by  DeRusso  (reference  6,  Chapters  1 through  4).  We  briefly  introduce 
these  topics  to  define  notation. 

It  can  be  shown  (reference  6,  p.  23)  that  a time-invariant  linear  differ- 
ential equation  has  an  output  of  y(t)  expressed  as  a function  of  the  Input  x(t) 
by  the  equation 


y(t) 


h(t-r)  x(t)dt 


h(r)  x(t-r)dt 


(2-56) 


where  h(-r)  Is  the  Impulse  response  of  the  system.  Such  a system  Is  depicted 
dlagrammatlcally  In  figure  12. 
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X(t)  w 
X(w) 


Figure  12.  A Linear  Dynamic  System. 

The  Fourier  transform  X(w)  of  a signal  x(t)  Is  a complex  function  representing 
the  magnitude  and  phase  of  the  frequency  components  of  x(t).  The  transform  Is 
defined  as 


F[xCt)]  i X(w) 


x(t)exp(-jwt)dt* 


(2-57) 


An  Inverse  transformation  also  exists  by  which  x(t)  my  be  obtained  from  X(w). 


x(t) 


X(w)  exp(jwt)dw 


(2-58) 


Observe  that  the  transform  Is  a linear  operator. 
EXAMPLE  2-4 


We  explicitly  calculate  the  Fourier  transform  of  x(t)  a A sin  (2irft)  to 
show  that  the  transform  represents  the  frequency  or  spectral  content  of  the  cor- 
responding time  signal.  Strictly  speaking,  the  Fourier  transform  of  this  par- 
ticular x(t)  does  not  exist  since  x(t)  Is  not  absolutely  Integrable,  l.e., 


|As1n(2irft)  | dt  * » 


*The  quantity  j Is  the  square  root  of  -1 , j = 
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Consequently  we  cannot  apply  equation  (2-57)  directly  to  evaluate  X(w).  Ob- 
serve, however,  that  assuming 


x(w)  = jirA[<$(w+2irf)  - <S(w-2irf)] 


where  6 Is  the  Dirace  delta  function,  we  can  verify  that  X(w)  Is  Indeed  the 
transform  of  x(t).  We  evaluate  the  Inverse  transform  by  use  of  the  delta  func- 
tion sifting  property  (reference  6)  to  obtain: 


X(w)  exp(jwt)  dw  = j j[exp(-j2irft)  - exp  (j2irft)l 


* A sin  (2irf t)  3 x(t) 


This  transform  pair  Is  shown  In  figure  13.  Observe  that  as  the  frequency  of 
the  sinusoid  Increases,  the  transform  delta  functions  move  farther  from  the 
origin.-* 


|X(»)| 


LI 

0 2vf  w 


Figure  13.  Fourier  Transform  Pairs. 

The  linear  dynamical  system  Input-output  relationship  given  In  equation 
(2-56)  can  be  Fourier  transformed  to  yield  an  equivalent  transform  description 
of  the  system.  Performing  the  transform,  we  have 


Y(w) 


■/: 


y(t)  exp  (-jwt)  dt 


■/;[/; 


h(t-r)  x (t)  dr 


exp  (-jwt)  dt 
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We  Interchange  the  order  of  Integration  and  multiply  the  Integrand  by  the 
quantity 

exp  (jwt)  exp  (-jwt)  * 1 

and  obtain 


h(t-x)  exp  [-jw(t-r)]  dtjexp  (-jw-r)  x(t)  dr 

= H(w)  X(w)  (2-59) 

Equation  (2-59)  shows  the  simple  relationship  that  the  Fourier  transform  of  the 
output  of  a linear  system  Is  the  product  of  the  Fourier  transforms  of  the  Input 
signal  and  the  system  Impulse  response  function.  A similar  Input-output  re- 
lationship is  now  shown  for  the  system  driven  by  stochastic  forcing  functions. 
Since  Fourier  transform  techniques  are  not  applicable  to  stochastic  process 
(they  are  not  absolutely  Integrable),  a few  preliminary  concepts  must  be  de- 
fined prior  to  the  formulation  of  a stochastic  process  transfer  function. 

Signal  frequency  distributions  for  stochastic  processes  are  derived  by  the 
following  heuristic  argument.  Given  any  stochastic  process  x(t)  we  may  filter 
that  process  with  a narrow  band  filter  having  center  frequency  f and  bandwidth 
Af.  The  transform  of  an  Ideal  bandpass  filter  is  plotted  In  figure  14.  The 


1 


Figure  14.  An  Ideal  Bandpass  Filter 
Transform  Characteristic. 
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output  signal  of  this  bandpass  filter,  just  as  In  the  analogous  case  In  which 
a deterministic  signal  Is  filtered,  represents  the  Input  signal  spectral  con- 
tent within  the  passband  frequencies.  The  power  of  the  filter  output  repre- 
sents that  portion  of  the  Input  signal  power  within  the  passband  of  the  filter. 
We  define  signal  power  P as  the  time  average  signal  squared  for  both  deter- 
ministic and  stochastic  signals.  This  time-average  for  stochastic  signals  Is 
the  autocorrelation  function  evaluation  for  zero  argument 


P = 


x2(t)  dt 


= mxx(0)  (for  stochastic  processes) 


(2-60) 


The  power  output  of  the  filter  can  also  be  measured  by  time-averaging  the  signal 
squared.  The  complete  analysis  procedure  of  bandpass  filtering,  squaring,  and 
time-averaging  Is  depicted  In  figure  15a.  We  denote  the  average  power  density 


fJV(t)] 


Figure  15a.  Measuring  Narrow  Band  Signal  Power. 
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Figure  15b.  Measuring  Narrow  Band  Cross-Spectral  Power. 


Cross-Spectral  Densities  (CSD),  * (f),  between  two  signals  x(t)  and  y(t)  are 
analogously  defined  but  must  be  generalized  in  order  to  characterize  the  Impor- 
tant phase  relationships  between  two  distinct  signals  that  do  not  arise  In  con- 
siderations of  a single  signal.  For  PSD  measurements,  a signal  x(t)  Is  passed 
through  a narrow-band  bandpass  filter  which  gives  a sinusoidal  output  signal 
y(t)  having  an  amplitude  proportional  to  the  square-root  of  the  signal  power 
within  the  filter  pass-band.  The  power  Is  determined  by  the  time-averaged  fil- 
ter output  squared.  The  cross -spectral  power  between  signals  x(t)  and  y(t)  Is 
calculated  by  multiplying  the  outputs  of  two  narrow  band-filters  excited  by 
, x(t)  and  y(t),  respectively;  then  the  time  average  product  Is  computed  as  shown 

• f 

\ In  figure  15b.  Note  that  this  procedure  gives  the  proper  PSD  In  the  case  that 

y(t)  equals  x(t).  For  more  general  signals  of  x(t)  and  y(t),  it  represents  the 
In-phase  signal  power  common  to  signals  x(t)  and  y(t).  An  out-of-phase  compo- 
nent, the  quadspectrum,  must  also  be  computed  to  complete  the  two  signal  common 
characterization.  The  requirement  for  calculation  of  the  quadspectrum  follows 
from  consideration  of  the  case  In  which  y (t)  Is  a delayed  or  phase-shifted  ver- 
sion of  x(t)  such  that  the  phase  shift  y(t)  with  respect  to  x(t)  Is  precisely 
90°.  In  this  case  the  averaged  Input  Is  the  product  of  orthogonal  signals  sine 
and  cosine  and  the  output  will  be  zero.  By  shifting  one  of  the  Input  signals 
90°  at  the  pass  band  center  frequency,  fQ,  the  roles  of  the  In-phase  cospectrum 
and  the  out-of-phase  quadspectrum  are  reversed.  Thus  by  Inserting  a 90°  phase 
shifter  In  the  y(t)  signal  path,  the  quadspectrum  between  x and  y Is  determined 
by  procedures  otherwise  Identical  to  those  employed  to  determine  the  cospectrum. 
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Denoting  the  cospectrum  power  by  P we  have  that 

c 


Pc(x,y,f0,Af)  * Rzw(o)/Af 


(2-63a) 


and  analogously  the  quadspectrum  power  Pq  Is 


pq(x,y,f0,Af)  - Rw(x)/Af 


(2-63b) 


where  t represents  the  time  shift  required  to  change  the  relative  phase  of  x 
with  respect  to  y exactly  90°  at  the  passband  center  frequency. 


27rfQT  = tt/2 


T'^ 


Rzw(t)  rePresents  cross  correlation  function  of  signals  z and  w obtained  by 

bandpass  filtering  signals  x and  y respectively.  The  cross-spectral  density 

» (f)  Is  defined  In  terms  of  the  co-  and  quad-spectrum  as 
xy 


*xy(f)  [Pc(x’y’Vf)  - jPq(x*y*f0*Af)] 

Af-t-0 


(2-64) 


Unlike  the  power  calculations  for  PSDs  which  must  always  yield  a non-negative 

result,  P and  P may  assume  any  real  value.  Also  note  that.  In  general,  the 
c q 

cross-spectral  density  Is  a complex  function  of  frequency.  These  generaliza- 
tions over  the  PSD  are  required  to  account  for  the  phase  relationships  betwen 
two  signals. 

PSDs  and  CDS  are  shown  to  be  directly  related  to  Fourier  transforms  of  the 
auto  and  cross-covariance  functions.  By  direct  calculation 
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Rxy(r)  exp  (-jwx)  dx 


* 11m 
T-*« 


T/Z 

VJ* 

T/2 


x(t+x)  y(t)  dt  exp(-jwx)  dt 


(2-65) 


For  mathematical  convenience  we  define  a truncated  function  Xy(t) 


xT(t) 


x(t)  |t|  <T/2 

0 otherwise 


(2-66) 


and  defining  y^(t)  similarly,  substitute  these  relations  Into  equation  (2-65) 
(and  Ignoring  negligible  end  effects  we  have) 


F(w) 


ilm  1/T 
T— 


xT(t+x)  yT(t)dt  exp(-jwx)dx 


(2-67) 


Multiplying  the  Integrand  by  exp  (-jwt)  exp  (+jw  ) = 1 and  rearranging  gives 


F(W) 


s '"DL 


Xj( t+x )exp[- jw ( t+x ) ]dt 
exp  (jwx)dx 


X-r(w)  Yt(-w) 

* an  -! 

1 


(2-68) 


Recall  from  equation  (2-58)  and  Example  2-4  that  the  Fourier  transform  represents 
a signals  spectral  content.  From  Parseval's  identity  for  Fourier  transform  pairs, 
namely. 
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x(t)  y(t)dt 


X(w)  y(-w)dw 


(2-69) 


we  see  that  the  Fourier  transform  magnitude  squared  represents  precisely  the  sig 
nal's  energy  distribution  In  frequency.  Specifically  we  have 


but 


x2(t)dt 


X(w)  X(-w)dw 


X(w)X(-w)  = |X(w) |2 


Consequently, 


(2-70) 


1 

27 


-2ir  (f0-Af/2) 


“2ir  (fft+Af/2) 


|X(w) [2dw  + 


21r(fQ+Af/2) 

|X(w) |2dw 

(f0-Af/2) 


represents  a signal  energy  In  a f frequency  band  about  center  frequency  fQ. 
These  concepts  for  signal  energy  properties  are  also  applicable  to  signal  power 
considerations  because  energy  normalized  by  time  Is  power.  Consequently,* 


*A  subtlety  In  the  concept  of  frequency  distribution  of  power  has  arisen  upon 
which  we  briefly  elaborate  here.  The  bandpass  conceptualization  of  power  In 
equation  (2-61)  and  equation  (2-62)  Is  valid  for  positive  frequencies  only. 
We  cannot  build  a filter  having  only  negative  bandpass  frequencies.  Any 
realizable  filter  passes  equally  frequencies  f2  = - fi.  Thus  we  define  PSD 
and  CSD  for  positive  frequencies  only.  However,  Fourier  transforms  define 
spectral  content  for  both  positive  and  negative  frequencies,  hence 


V,(f)df 


and  equations  (2-71)  and  (2-72)  follow. 


48 


AFWL-TR-76-193 


Xt(w)  Yt(-w) 
4„>/2h)  s 2 Zim  - =? 

* TV« 

3 2 F [Rxy(r)]  i 2 Gxy(w) 


(2-71) 


as  a special  case 


# % |xT(w)|2 

®yy(w/2n)  = 2 lim  — L 

xx  T-*»  1 


= 2 F [R(t)]  * 2 G (w) 

xx  xx  w > 0 


(2-72) 


The  characterization  tools  that  have  now  been  developed  allow  us  to  calculate 
Input-output  signal  relationships  for  linear  dynamic  systems.  Equation  (2-56)  . 

repeated  below  describes  a linear  system  output  signal  y(t)  in  terms 

of  the  system  impulse  response  h and  the  system  input  x.  This  equation 

is  equally  valid  for  stochastic  process  excitation  of  the  system. 


h(rj)  x(t-Tj)  dTX 


(2-73) 


Taking  expectations  of  equation  (2-73)  gives  the  system  output  mean  as  a function 
of  the  input  mean.  For  stationary  processes  the  output  mean  is  simply  a scalar 
multiple  of  the  input  mean. 


uy  «£[y(t)] 


h(xi)ux  df! 


h(x i ) dt! 


* ux  H(w) 


w = 0 
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Hence,  the  average  value  of  the  system  output  Is  simply  the  average  value  of  the 
system  Input  multiplied  by  the  dc  gain  H(o)  of  the  linear  system. 

Higher  order  moments  of  the  system  output  and  joint  moments  of  the  system 
output  and  Input  are  similarly  obtained  as  functions  of  the  linear  system  Im- 
pulse response  function  h(t)  and  moments  of  the  Input  signal.  For  example,  the 
cross-correlation  of  the  system  input  and  output  signals  Is  obtained  as  follows. 
Multiplying  both  sides  of  the  equation  (2-73J  by  x(t  + t)  and  taking  expecta- 
tions yields: 


Vr) 


h(Ti)  Rxx(t+ti)  dxi 


(2-74) 


This  equation  shows  that  the  cross-correlation  of  the  system  input  and  output 
signals  Is  given  as  the  convolution  of  the  system  Input  signal  auto-correlation 
function  with  the  system  Impulse  response.  Taking  the  Fourier  transform  of  both 
sides  of  equation  (2-74)  gives: 


4>xy (W/2lr)  = H("W)  *xx(W/2lr) 


(2-75) 


Similarly,  multiplying  both  sides  of  equation  (2-73)  by  y(t  + t)  and  taking  ex- 
ceptions gives: 


V°  “ 


hUl)  RyX(T+Tl)dT! 


(2-76) 


and  transforming: 


*yy(w/2ir)  = H(-w)  *yx(w/2u) 


* H(w)  *xy(w/2ir) 


(2-77) 


By  the  definition  of  correlation  functions  and  their  transforms  we  have  the  fol- 
lowing properties: 


50 


AFWL-TR-76-193 


V°  - Ryx(-x) 


*xy(w/2ir)  = *yx(w/2ir)  = $xy(-w/2w) 


(2-78) 


Thus  substituting  for  * In  equation  (2-77)  from  equation  (2-75)  by  using  equa- 
tlon  (2-78)  gives  the  desired  result. 


$yy(w/2ir)  = H(-w)  *xy(w/2ir) 


- H(-w)[H(-w)*xx(w/2ir)]* 


= | H (w)  | 2 *xx(w/2ir) 


(2-79) 


since  $xx  is  real. 


We  also  write 


«yy(t)  = |H(2irf)|2  $xx(f) 


(2-80) 


’’his  equation  shows  that  a linear  dynamical  system  output  signal  PSD  Is  precisely 
the  Input  PSD  'spectrally  shaped1  by  the  frequency  domain  transfer  function  of  the 
system.  We  illustrate  the  response  of  a linear  dynamical  system  excited  by  a 
stochastic  process  with  an  example. 

EXAMPLE  2-5 

A simplified  linear  dynamic  system  is  considered.  The  system  Is  modeled  by  the 
first  order  linear  time-invariant  differential  equation 


ay  + ax 


The  quanlty  y represents  the  system  output  while  x represents  a stochastic  pro- 
cess excitation  of  the  system.  The  differential  equation  specified  above  repre- 
sents a lowpass  filter.  Spectral  components  of  x having  frequencies  f less  than 
the  filter  cutoff  frequency  fc  = a/2ir  are  passed  unattenuated  to  the  output  y. 
Spectral  components  of  x at  frequencies  greater  than  fc  are  attenuated.  The 
linear  system  Impulse  response  Is 
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h(t)  = a exp(-at)  t ^ o 

and  Its  Fourier  transform,  also  known  as  the  system  transfer  function,  is 

H(w)  = (1  + jw/a)-1 

We  specify  that  the  stochastic  excitation  is  a zero  mean  process  with  an  auto- 
correlation function 


Rxx(0  = ^(r) 

This  stochastic  process  is  of  particular  importance  in  analysis.  Since  Its 
autocorrelation  function  Is  zero  for  all  t not  equal  to  zero,  we  conclude  that 
knowledge  of  the  value  of  x at  time  t is  uncorrelated  with  its  value  at  any 
other  time.  Hence  we  cannot  predict  future  values  of  x(t)  from  past  observa- 
tions. Stochastic  processes  having  autocorrelation  functions  of  the  form  q<$(t) 
are  called  white  noise.  This  terminology  follows  from  its  spectral  density, 
which  is  constant  for  all  frequencies,  and  the  analogy  with  white  light,  which 
Is  composed  of  light  of  all  colors  or  frequencies,  the  excitation  spectral 
density  Is  derived  from  equation  (2-72)  and  Rxx(0  as 

= 2q  ; wSso 


Using  equations  (2-75)  and  (2-80)  we  calculate  that 

= H(-w)  ^(w/Sit) 

= 2q  (1  - jw/a)-1 
and 


*yy(w/2ir)  = lH(W)|2  $XX(W/2,r) 

= 2q  [1  + (w/a)2]"1 


the  correlation  functions  R (t)  and  R (t)  may  be  calculated  in  two  ways: 

xy  yy 
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- 

i 


(1)  their  Fourier  transforms  may  be  deduced  from  * and  * and  Inverted 


or 


xy 


yy 


(2)  they  can  be  calculated  directly  from  h(t)  and  Rxx(t)  by  use  of  equa- 
tions (2-74)  and  (2-76).  We  shall  calculate  Rxytx)  and  Ryy(x)  by  both 
methods  to  Illustrate  methodology  and  to  demonstrate  the  equivalence  of 
the  alternative  approaches.  Since 


*xy(w/2ir)  = 2q[l-jw/a]_1 


we  conclude  that 


P[Rxy  (t)]  = qD-jw/a]-1 


which  we  Inverse  transform  (consult  reference  6 or  any  convenient  table  of 
Fourier  transform  pairs)  to  obtain 


Rxy(t)  = aq  exp  (a-r)  x -£0 

Observe  that  Rxy(x)  = 0 for  t > 0.  This  follows  directly  from  the  causalty  of 
the  system  output  which  can  be  correlated  only  with  past  values  of  the  Input. 

Similarly 

F CR,y(x) J = q[i+(w/a)2]-1 

=■•  q/2[(l+jw/a)_1  + (1-jw/a)*1] 
and 

Ry/(t)  = q|  exp(-a |t | ) 

By  the  alternative  method  we  calculate 
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(DIRAC  DELTA)  4 


R„(r) 


Rxx(t)*qS(T) 


Figure  16.  Linear  System  Input  and  Output  Autocorrelation  Functions. 
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The  corresponding  power  spectral  densities  are  plotted  in  figure  17.  The  param- 
eter 1/a  represents  the  correlation  time  of  the  output  stochastic  process  and 


Figure  17.  Linear  System  Input  and  Output  Power  Spectral  Densities. 

Is  Identical  to  the  time  constant  of  the  linear  system  when  excited  by  white 
noise.  *-* 

The  Important  relationships  developed  in  this  section  are  summarized  in 
table  1.  Also  Included  In  the  table  are  other  important  equalities  which  the 
reader  Is  Invited  to  validate.  These  equations  play  a fundamental  role  through- 
out the  remainder  of  this  guide. 


Table  1 

CORRELATION  FUNCTIONS  AND  SPECTRAL  DENSITY  RELATIONSHIPS 


Definitions 

R«(,)  " £ 5Tj 

/ x(t+x)  x(t)  dt 

' -T 

Autocorrelation 

VT)  = £ 5TJ 

j"  x(t+r)  y(t)  dt 
-T 

Crosscorrelation 

*xx(w/2ir)  s 2 FCRxx 

(t)1  w^O 

PSD 

*xy(w/2^)  = 2 FCR^Ct)]  w^O  CSD 
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Spectral  $yv(f)  = H(-2irf)*  (f) 

Density  y 

Relationships 

« H*(2irf)*xx(f) 

*yx(f)  - H(2irf)*xx(f) 

= H*(-2irf)*xx(f) 

*yy(f)  = |H(2irf ) 1 2*xx(f ) 
♦Denotes  complex  conjugate. 


9 TRANSFER  FUNCTION  ESTIMATES  AND  COHERENCE 

In  the  previous  section  the  relationships  between  stochastic  process  prop- 
erties of  the  Input  and  output  signals  of  a linear  dynamic  system  are  derived, 
given  the  properties  of  the  linear  system.  Often  in  engineering  practice  we 
suspect  that  two  observed  stochastic  signals  are  related  by  some  linear  system 
which  we  have  not  characterized.  In  this  section  we  are  concerned  with  the 
questions:  Are  two  observed  signals  related  by  a linear  dynamic  system?  And 
If  so:  What  Is  the  transfer  function  of  that  system? 

Given  any  two  time  functions  x(t)  and  y(t),  we  could  compute  their  power 
and  cross-spectral  densities,  and  then  by  using  equations  of  table  1 solve  for 
an  estimate  H(w)  of  the  unknown  transfer  function.  This  procedure  gives 


H(w)  * 4yx(w/2ir)/ «xx(w/2ir)  w^O 
and 

H(w)  * H*(-w)  w < 0 (2-81) 

If  only  the  magnitude  of  the  estimated  transfer  function  were  required,  an 
alternative  approach  is 
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| H(w)  | 2 ■ *yy(w/2ir)/  $xx(w/2ir) 


(2-82) 


or  equivalently 


H(w)  = .J*  (w/2-ir)  (w/2ir) 


(2-83) 


As  subsequently  derived,  equations  (2-81)  and  (2-82)  give  the  best  possible 

estimates  (In  the  mean-square-error  sense)  available  from  the  data  y(t)  and  x(t). 

However,  additional  calculations  must  be  made  to  ascertain  whether  or  not  the 

signals  x(t)  and  y(t)  are  indeed  related  by  a linear  system.  These  calculations 

are  required  since  the  quantities  * and  $ are  strictly  positive  and  hence 

yy  xx 

nonzero  transfer  functions  estimates  H(w)  are  obtained  for  any  two  signals  x(t) 
and  y(t)  regardless  of  their  connection  to  a single  linear  system.  The  signals 
need  not  even  be  measured  simultaneously  to  give  non-zero  transfer  function 
estimates.  We  shall  define  and  compute  a coherence  function  which  plays  a role 
analogous  to  that  of  the  correlation  coefficient  introduced  in  section  II. 5 for 
random  variables  to  characterize  the  extent  of  linear  correlation  between  two 
time  signals. 

Recall  from  section  II. 5 that  given  data  sets  x^  and  y^  an  estimate  y^  can 
be  derived  which  Is  a linear  function  of  x^.  The  coefficients  a and  b of  the 
functional  relationship  given  In  equation  (2-84)  are  derived  to  minimize  the 
mean-square-error. 


yj  = axj+b 


(2-84) 


As  shown  In  section  II. 5 and  Example  2-3,  the  coefficients  a and  b can  be  calcu- 
lated from  the  covariance  and  the  variances  of  the  random  variables  x and  y or 
they  can  be  estimated  from  the  data  sets  themselves.  An  analogous  procedure  Is 
derived  from  estimating  the  linear  dependence  of  two  stochastic  processes  x(t) 
and  y(t). 

Given  the  stochastic  signals  x(t)  and  y(t),  we  assume  that  they  are  related 
by  a linear  system.  Logically,  an  estimate  y of  y can  be  written  as  a convolu- 
tion Integral  like  that  of  equation  (2-56).  Explicitly  we  form  the  estimate 
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. 

I 

I 


h(x)  x(t-r)  dr 


(2-85) 


Since  only  the  functions  x(t)  and  y(t)  are  given,  no  prior  knowledge  of  the  sys- 
tem Impulse  response  function  h(x)  or  Its  Fourier  transform  H(w)  Is  available. 
These  functions  must  be  estimated  from  the  given  data.  Some  estimation  error 
criterion  must  be  formulated  and  the  transfer  functions  chosen  to  minimize  the 
error  functional.  The  mean-square  error  functional  leads  most  simply  to  calcu- 
lable results.  The  estimation  error  is  defined  as 

e(t)  = y(t)  - y(t)  (2-86) 

The  autocovariance  of  the  error  Is  readily  calculated  as 


Ree(r)  «£[e(t+x)  e(t)] 


h ■/.: 


h(x! ) x(t+T-Tx)  dxX 


1 


h(x2)  x(t-x2)  dx2 


- 

> 


•r 


Ryy(T)  - I h(xj)  [^(-Ti-x)  + Rxy(T“Tl)]  dTl 


h(xi)  h(x2)  rxx(xi-x2  -x)  dxj  dx2 


(2-87) 


The  mean-square  error  can  be  minimized  by  proper  choice  of  the  system  impulse 

response  h(t).  We  define  F as  those  components  of  Rge  (0)  dependent  upon  h. 

We  minimize  RQO(0)  by  equivalently  minimizing  F. 
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•f 


h(^l)  RxyHl)  dTl 


not) 
00 


h(xx)  h(x2)  Rxx(ti-t2)  dt!  dx2 


■/:  ■""[/; 


h(t2)Rxx(T1  -t2)  dx2  - 2 Rxy(-Tl)  | dTl 


] 


(2-88) 


A variational  method  yields  conditions  which  the  minimizing  function  h(x)  must 
satisfy.  Let  h°(x)  denote  the  minimizing  solution  and  form 


h(x)  = h°(x)  + eh^x) 


where  h!(x)  Is  an  arbitrary  function  and  e is  a scalar  parameter.  Substituting 
for  h(x)  Into  equation  (2-88)  and  simplifying  gives 


F = F°  + eF1  + e2  F2 


(2-89) 


where 


F° 


h®(x2)  Rxx(t:-t2)  dx2  - 2 Rxv(-xi)  dxj 


xy' 


■] 


(2-90) 
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l 


F 


hl  (Tl)  RXy(-Ti)  dTi 


h°  (t2)  +h°(xi)  h1(r2)]  Rxx(t1  _x2)  dxj  dx2 


00 

h°(T2)  Rxx(t1  -*z)  dT2  - RXy(“T 1 ) 

00 


dxi 


(2-91) 


and 


F2 


h1 (t i ) h1 (t2)  rxx(ti  -t2)  dtj  dx2 


(2-92) 


It  can  be  shown  that  F2^0  for  arbitrary  h^x). 

The  F1  term  dominates  any  variation  of  F with  e for  e sufficiently  small.  Since 
e may  assume  both  positive  or  negative  values,  a necessary  condition  that  h°  (x) 
be  the  minimizing  solution  Is  the  requirement  that 


F1  = o for  arbitrary  hx(x) 
which  In  turn  requires  that 


V'T) 


h°  (T2 ) Rxx(t_  t2 ) dx2 


(2-93) 


(2-94) 


The  Integral  equation  (2-94)  Is  one  variant  of  the  so-called  Wlener-Hopf  equa- 
tion (reference  5,  p.  305).  Solution  for  a function  h°  (x2)  which  satisfies 
equation  (2-94)  Is  most  readily  obtained  by  Fourier  transforming  that  equation 
to  obtain  an  equation  for  H°(w)  which  can  be  Inverse  transformed  to  yield  a 
solution.  Effecting  the  transform 
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F[Rxy(-x)]  = H°(w)  F[Rxx(tx)] 


(2-95) 


denoting 


Gyx(w)  * F[Ryx  (t)]  = F [Rxy  (-x)] 


and 


Sxx<»)  *FCRxx<’H 


gives  the  solution 


H°(w)  = Gyx(w)/Gxx(w) 


and 


h°(r)  = F-l[H°(w)] 


(2-96) 


(2-97) 


(2-98) 


(2-99) 


In  general  the  Impulse  response  function  obtained  by  this  technique  is  nonzero 

for  t < 0.  Thus,  no  physical  system  could  be  constructed  having  this  Impulse 

function  since  all  realizable  systems  are  non-anticipating  (they  cannot  begin  to 

respond  to  an  Input  that  has  not  yet  been  applied).  A spectral  factorization 

technique  developed  by  Norbert  Wiener  (reference  7)  does  generate  a solution  of 

equation  (2-94)  which  Is  Implementable  by  casual  systems.  This  approach  is  most 

readily  used  In  the  case  that  the  functions  G (w)  and  G (w)  are  expressed  as 

yx  xx 

rational  polynomials  In  w.  Since  such  polynominal  expressions  are  not  avail- 
able from  our  analysis  we  discuss  his  technique  no  further. 

Summarizing  these  results,  If  Gyx(w)  = F [Rxy  (-x)]  and  Gxx(w)  = F [Rxx  (r )], 
then  the  best  estimate  y(t)  of  y obtainable  from  the  data  x(t)  is 


h(r)  x(t-t)  dt 


(2-100) 
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or  equivalently 

dw 

(2-107) 

This  latter  expression  is  particularly  useful  since  it  allows  a definition  of 
normalized  error  on  a spectral  basis.  We  define  the  coherence  function  of  the 
estimate  as 


^ -o 


V> 


1 |S»»12 


xy 


(f) 


lGYX(w)l2  l«vx(0!2 

G (w)  G (wT  " $ (f ) 4 (TJ 
xx'  yy'  ' xxv  ' yyv  7 


(2-103) 


Clearly 


0 — Yxx(f ) — ! 


(2-109) 


The  coherence  function  plays  a role  similar  to  the  correlation  coefficient.  It 

characterizes  the  degree  of  the  linear  relationship  between  two  signals  such 

that  as  YyV(f ) approaches  1 the  tv/o  functions  are  more  highly  correlated,  while 

yMf)  equal  to  zero  guarantees  no  linear  dependence, 
xy 

EXAMPLE  2-6 


Consider  a linear  dynamic  system  with  input  x(t)  and  response  y(t).  Power 
spectral  densities  of  the  input  and  output  signals  have  been  measured  as  well  as 
the  Input-output  cross  spectral  density.  An  estimate  of  tne  system  transfer 
function  is  to  be  obtained  from  these  quantities.  The  coherence  function  is 
also  computed  and  the  results  interpreted.  The  measured  spectral  densities  are 


Sxx(w)  ■ 1 


V"1  ' l+(w?aj-  * -0' 


and 


Gyx^  " 1 + j w/a 
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From  equation  (2-98)  the  transfer  function  estimate  is  computed  as 

Gvy(w)  i 
H(w)  = = Ti^57i’ 


The  coherence  function  Is  computed  by  substitution  for  G (w),  G (w),  and 
Gxx(w)  in  equation  (2-108).  yx  yy 

Y2  (w/2»)  - I1  * Wal~2  - Jl  + .01  [1  + (w/a)2] 

Xy  T-rW  + -0’  ' 


Observe  that  at  low  frequencies  where  w < a we  have  *r2  » 1.  At  these  frequen- 
cies the  estimate  H(w)  of  the  transfer  function  is  very  good.  At  frequencies 
w > a,  yz  Is  less  than  one  and  it  approaches  zero  as  w becomes  large.  The 
estimate  is  not  good  at  high  frequencies.  A block  diagram  description  of  a 
system  having  the  properties  given  in  this  example  is  presented  in  figure  18. 


Figure  18.  Linear  System  Transfer  Function  Identification  for  Example  2-6. 


1 he  estimated  transfer  function  agrees  identically  with  the  actual  system  trans- 
fer function.  The  coherence  function  decreases  from  unity  at  higher  frequencies, 
however,  because  the  measurement  noise  n(t)  becomes  proportionately  large  with 
respect  to  the  linear  system  output.  At  very  high  frequencies,  the  signal  y(t) 

Is  dominated  by  n(t)  and,  consequently,  y(t)  is  not  accurately  modeled  by  y(t) 
as  obtained  from  equations  (2-100)  and  (2-101).  The  coherence  function  devia- 
tion from  unity  Indicates  this  fact.#-* 

Additional  properties  of  the  coherence  function  and  its  role  in  transfer 
function  estimation  are  developed  in  section  IV.  We  caution  the  reader  that 
these  techniques  are  applicable  to  stochastically  excited  systems.  Transfer 
function  estimates  for  linear  systems  with  deterministic  input  signals  are  ob- 
tained by  other  techniques. 
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Other  useful  formats  for  the  presentation  of  PSD  information  have  been  de- 
vised. These  formats,  all  entailing  integrals  of  PSD!s  over  a given  frequency 
band,  are  used  most  commonly  for  characterizing  the  performance  of  stochasti- 
cally excited  systems.  Definitions  for  cumulative  power  (CUM  PWR),wide  band 
RMS,  a.  Inband  RMS,  and  BACKWARD  SUM  are  presented  here  for  completeness.  We 
comment  briefly  upon  the  Interpretation  of  each  term. 


CUM  PWR  (fH)  « 


(2-110) 


Cumulative  power  represents  a signal's  mean  square  value  in  frequency  components 
from  dc  to  a maximum  frequency  f^.  Wideband  RMS  represents  the  root-mean-square 
signal  value  contained  In  signal  frequency  components  from  dc  to  a maximum  fre- 
quency. RMS  Is  the  level  of  a dc  signal  containing  the  same  root-mean-square 
value  as  the  time-varying  signal  being  analyzed. 


Wide  Band  RMS  = 


CT  = 


*xx(f)  df 


(2-111) 


Clearly, 


o 


^COM  PWR  (-) 


Inband  RMS  is  simply  a generalization  of  wideband  RMS.  The  lower  and  upper  fre- 
quency limit  are  specified.  Inband  RMS  represents  the  signal  power  within  the 
specified  frequency  band. 


In  Band  RMS  (fL>  fH)  = 


Vf> " 


(2-112) 


Finally,  BACKWARD  SUM  (f^)  represents  the  signal  power  attributable  to  frequency 
components  above  the  specified  frequency  f^. 
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BACKWARD  SUM  (f 


■/; 


$xx(f)  df 


= a2  - CUM  PWR  (fL) 


10  SUMMARY 

The  key  concept  developed  in  this  chapter  Is  that  of  a stochastic  process 
representing  signals  and  phenomena  which  cannot  be  accurately  predicted  and  are 
of  a random  or  statistical  nature.  These  signals  are  characterized  in  terms  of 
their  mean,  covariance  and  spectral  density  properties.  Given  these  properties 
of  a linear  system  input  forcing  function,  and  the  system  transfer  function,  one 
can  calculate  the  corresponding  properties  of  the  system  response.  These  calcu- 
lations are  important  when  calculating  the  response  of  an  hypothesized  system 
(no  hardware  exists)  to  a known  stochastic  environment  and  for  the  performance 
evaluation  of  such  systems.  The  transfer  function  estimation  techniques  pre- 
sented in  section  II. 9 allow  calculation  of  the  linear  dependence  of  two  sig- 
nals. These  techniques  can  establish  cause  and  effect  relationship  between 
seemingly  unrelated  signals. 
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SECTION  III 

COMPUTATIONAL  ASPECTS  OF  ANALYSIS 


1 INTRODUCTION 

The  primary  emphasis  of  this  section  is  the  development  of  specific  data 
analysis  procedures  based  upon  the  theoretical  results  presented  in  the  last 
section.  In  developing  practical  data  analysis  techniques,  we  shall  prefer  (or 
in  some  cases  be  forced)  to  make  simplifying  approximations.  These  approxima- 
tions shall  be  examined  in  detail  to  determine  their  Impact  upon  the  interpre- 
tation of  our  results.  The  basic  results  of  this  section  deal  with  development 
of  algorithms  suitable  for  use  with  finite  duration,  sample-data  sequences. 

These  numerical  procedures  are  suitable  for  coding  in  a higher  level  language 
and  execution  in  a minicomputer  or  a larger  computer  system.  The  Discrete 
Fourier  Transform  (DFT)  is  introduced  as  a sampled-data  equivalent  of  the 
Fourier  transform.  The  Fast  Fourier  Transform  as  a computationally  efficient 
algorithm  for  calculating  the  DFT  of  a numerical  sequence  is  explained.  Equa- 
tions for  computing  power  and  cross-power  densities  and  the  coherence  function 
are  presented  at  the  conclusion  of  this  section. 

2 FINITE  DURATION  DATA  INTERVALS 

The  power  spectral  density  of  a stochastic  process  x(t)  and  the  cross- 
spectral  density  of  two  processes  x(t)  and  y(t)  have  been  shown  to  be  Important 
characterizations  of  stochastic  processes.  These  properties  are  sufficient  for 
estimation  of  the  linear  dependence  of  two  processes  and  for  calculating  the 
output  power  spectrum  of  a dynamic  linear  system  excited  by  a stochastic  process. 
In  section  II  the  power  spectrum  is  shown  to  be  the  Fourier  transform  of  the 
process  autocorrelation  function,  equation  (2-71).  Thus  the  process  PSD  is 
readily  calculated  once  the  autocorrelation  function  has  been  determined. 
Theoretically,  the  autocorrelation  function  is  determined  by  the  second-order 
density  function,  equation  (2-44),  or,  for  ergodlc  processes,  by  the  time- 
average  correlation  computed  from  an  observed  realization  of  the  process,  equa- 
tion  (2-51).  Second-order  density  functions  are  generally  not  known  a priori. 
Estimation  of  a process  autocorrelation  function  is  generally  much  simpler  than 
estimation  of  the  second-order  density  function  and  calculation  of  the  auto- 
correlation function  from  it.  In  either  approach,  we  are  practically  limited 
by  a finite  duration  interval  of  data  from  which  the  estimates  can  be  computed. 

We  explicitly  examine  the  effects  of  the  finite  duration  data  interval  in  the 
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calculation  of  a power  spectral  density  as  the  Fourier  transform  of  a time  aver- 
age  autocorrelation  function  estimate. 

The  finite  duration  Interval  has  two  Immediate  consequences  in  the  Inter- 
pretation of  our  results.  Limiting  the  data  interval  to  a finite  duration 
makes  the  autocorrelation  function  a random  variable  with  a nonzero  covariance. 
The  finite  Interval  also  distorts  the  estimated  spectrum  of  Rxx(f). 

The  randomness  Introduced  by  using  limited  data  is  best  described  by  the 
following  observations.  Define 


Xj( t ) Xj(t  +t)  dt 


(3-1) 


where 


xT(t)  = 


x(t)  , o — t ^ T 
o,  otherwise 


(3-2) 


Note  that  RXyXj(f)  is  certainly  a random  variable  since  it  is  a function  of  the 
stochastic  processes  x(t)  and  x(t+t).  Under  a general  set  of  assumptions  (ref- 
erence 4,  section  9)  the  covariance  of  RXjX^(t)  approaches  zero  as  T approaches 
Infinity.  For  engineering  purposes  we  can  treat  the  estimated  correlation  func- 
tion just  like  any  other  deterministic  function  of  t.  The  finite  duration  inter- 
val, however  restricts  ‘I  to  finite  values.  Hence,  the  covariance  of  the  random 
variable  RXyX^.(T)  remains  non-zero.  The  finite  covariance  of  the  autocorrela- 
tion function  estimate  has  the  following  practical  interpretation.  If  we  were 
to  compute  RXjX^(t)  from  different  data  intervals  (or  from  different  members  of 
the  time  ensemble)  of  the  same  ergodic  stochastic  process,  we  would  calculate 
different  functions  RXjXT(t)  for  each.  The  mean  or  average  value  of  the  differ- 
ent calculations  is  Rxx(t);  the  various  values  are  dispersed  about  the  mean  with 
a covariance  which  (in  principle,  at  least)  is  calculable  from  the  density  func- 
tions of  the  underlying  probability  space.  A PSD  calculated  by  Fourier  trans- 
forming Rxtxt(t)  necessarily  has  a corresponding  covariance  or  uncertainty. 

Bounds  on  the  PSD  uncertainty  are  derived  in  section  III. 4. 

Calculation  of  a PSD  estimate  from  a finite  duration  data  interval  is  also 
claimed  to  have  spectral  distortion.  By  direct  calculation 
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g(t)  exp  (-jwt)  dt 


(3-8) 


Equation  (3-7)  Is  Interpreted  as  a complex  convolution  of  the  functions  X(w)  and 
G(w)  (compare  with  the  convolution  of  time  functions  given  In  equation  (2-56)). 

We  see  that  X^(w)  represents  a complex  weighted  average  of  the  values  X(wJ  In  a 
neighborhood  of  frequencies  about  w.  The  weighting  function  G(w)  is  simply  the 
Fourier  transform  of  the  data  window  function.  Generally  |G(w)|  has  largest  mag- 
nitude for  w « o and  Its  magnitude  approaches  zero  as  w becomes  large.  Typically 
a tradeoff  exists  between  the  main  lobe  width  and  the  rate  at  which  |G(w)|  rolls- 
off  with  Increasing  frequency.  A sharp  main  lobe  is  desired  so  that  Xy(w)  closely 
approximates  X(w).  A fast  rolloff  Is  required  to  prevent  components  X(wx)  from 
significantly  contributing  to  XT(w)  whenever  w^.  Examples  of  various  window 
functions  and  their  corresponding  transform  magnitudes  are  plotted  In  figure  19. 

The  spectrum  of  Xy(w)  Is  a distorted  representation  of  the  original  signal 
spectrum.  This  distortion  Is  most  readily  seen  *or  the  case  wherein  X(w)  has 
discrete  frequency  singularities,  as  for  example,  does  x(t)  = cos  wQt.  The 
actual  spectrum  |X(w)|  and  the  truncated  signal  spectrum  |Xy(w)|  are  plotted  In 
figure  20.  Notice  that  |Xy(w)|  has  non-zero  values  over  frequency  bands  where 
|X(w)|  Is  zero.  This  property  of  | Xy (w) | is  due  to  a leakage  phenomenon  associ- 
ated  with  all  finite  duration  window  functions.  The  primary  limitation  of  data 
windows  Is  that  they  smooth  the  measured  spectra,  introduce  leakage  components, 
and  limit  the  resolution  of  the  final  spectral  estimate.  These  effects  are  mini- 
mized by  a choice  of  T sufficiently  large  that  these  attendant  degradations  are 
inconsequential.  Generally  this  requirement  means  that  T > 2ir/Af  where  Af  Is 
any  frequency  difference  we  wish  to  resolve.  The  confidence  of  our  estimate 
will  also  require  larger  T for  better  estimates.  These  relationships  are  pre- 
sented In  Section  1 1 1. 4. 

The  cosine-taper  data  window  shown  in  figure  19d  Is  employed  in  the  AFWL 
SIGANAL  data  analysis  program.  This  window  function  provides  good  sidelobe  re- 
jection while  Introducing  minimum  distortion  to  the  original  data  signal.  A 
dominant  effect  of  the  distortion  is  to  reduce  the  variance  of  the  windowed  data 
with  respect  to  the  original  unwindowed  signal  (ref.  2,  p.  323).  The  windowed  data 
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Data  Windows  and  Their  Transform  Magnitudes 
fa)  Boxcar,  (b)  Triangular,  (c)  Hanning, 

(d)  Cosine  Taper. 


must  be  multiplied  by  a "correction  factor"  which  increases  the  output  data  vari- 
ation to  its  original,  unwindowed  value.  The  approximate  correction  factor  is 
given  by  the  ratio  of  the  window  function  'powers.'  A window  power  Is  simply  the 
area  under  the  curve  of  the  window  function  squared. 


CF 


^Boxcar  ^ 


dt 


1 

0.875 


'l 

(3-9) 


This  correction  factor  is  precisely  correct  for  random  data  as  can  be  seen  from 
the  following  calculations: 


while  for  windowed  data 
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(x(t)  - Px(t))2  g2(t)  dt 


€ 


g2  (t)  dt 


(3-11) 


The  ratio  of  these  variances  is  the  correction  factor  we  must  apply  to  approxi- 
mately boost  the  variance  of  the  windowed  data  to  match  the  original  signal  vari- 


g2(t)  dt 


(3-12) 


L_ 


An  alternant  method  that  can  be  used  to  calculate  the  approximate  correction  fac- 
tor Is  to  actually  calculate  the  sample  variances  of  data  both  before  and  after 
multiplication  by  the  window  function.  The  approximate  correction  factor  is  then 
simply  the  ratio  of  these  results.  Koenigsberg  (reference  8,  p.  80)  gives  an  ex- 
cellent discussion  of  this  procedure.  His  empirically  founded  observations  pro- 
vide useful  Interpretations  regarding  the  systematicness  of  the  data  wher  the 
correction  factors  as  calculated  above  differ  greatly  from  that  derived  from 
equation  (3-12). 

The  correction  factor  is  a power  boost  multiplier  for  signal  variances; 
hence.  It  Is  also  the  appropriate  correction  factor  for  PSD  calculations.  The 
complete  PSD  formulas  are  presented  In  section  III. 3. 

3 SAMPLED  DATA 

The  second  most  Important  aspect  of  practical  stochastic  data  analysis  pro- 
cedures deals  with  the  sampled-data  nature  of  signals  processed  by  numerical 
algorithms  on  a digital  computer.  The  equations  presented  thus  far  have  assumed 
that  values  of  the  time  function  x(t)  are  known  for  any  t within  the  analysis 
interval  [0,T].  Practically,  this  Is  the  case  only  when  the  signals  are  speci- 
fied analytically  in  a theoretical  deviation  or  when  the  data  analysis  equations 
are  Implemented  with  analog,  rather  than  digital,  circuits.  Prior  to  the  wide- 
spread use  of  digital  computers,  analog  bandpass  filter  techniques  were  employed 
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to  estimate  spectral  densities.  (Recall  the  PSD  discussion  In  section  II. 8.) 
Digital  numerical  procedures  make  feasible  and  practical  many  new  analysis  tech- 
niques. Consequently,  digitization  of  analog  data  signals  and  the  numerical 
evaluation  of  our  analysis  equations  requires  careful  examination  as  to  the  Im- 
pact of  these  operations  on  our  overall  results.  We  shall  find  that  data  digi- 
tization Introduces  distortions  of  the  estimated  PSD  just  as  did  the  finite 
duration  data  window. 

Data  digitization  or  sampling  typically  is  a process  whereby  a continuous 
signal  x(t)  Is  converted  to  a sequence  of  digital  words  represented  In  a binary 
number  system  common  to  modern  digital  computers.  Signal  values  corresponding 
to  the  digital  sequence  are  obtained  from  the  original  signal  at  fixed  time  Inter- 
vals aT.  The  sampling  frequency  or  sampling  rate  is  f = aT-1;  aT  is  the  sampling 
Interval  or  period.  A pictorial  representation  of  the  sampling  process  Is  pre- 
sented In  figure  21 . 


At  2 At  5 At  ioAt 


10  ' 


Figure  21.  The  Sampling  Process  for  Signal  Digitization. 

Assume  that  the  sample  points  X|  s x (1  AT)  for  integer  i are  known  rather 
than  the  values  of  the  function  x(t)  for  an  arbitrary  t.  We  wish  to  use  these 
sample  points  to  determine  Fourier  transform  spectral  properties  of  the  unsampled 
time  signal  x(t).  We  observe  that  one  possible  convention  for  reducing  the 
original  function  x(t)  to  a time  function  completely  characterized  by  the  values 
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x.|  Is  to  define  the  sampled-data  time  function  x*(t)  as  the  multiple  of  x(t)  with 
a time  sequence  of  Impulse  functions  spaced  aT  seconds  apart.*  The  sequence  of 
impulse  functions  Is  called  an  Impulse  train  and  is  denoted  6y(t). 


The  Fourier  transform  X*(w)  of  x*(t)  defined  by  equation  (2-57)  is 

I 

I 


♦ReceH  that  the  Impulse  function  (Dirac  delta  function)  Is  defined  by  the  follow- 
ing properties: 
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X*(w) 


*(t)  exp  (-jwt)  dt 


oo 

X*(w)  = x^xp  (-jw  aTI) 


(3-15) 


Certainly  we  could  omit  the  Intervening  steps  and  simply  define  a sample-data 
Fourier  transform  as 


oo 

X*(w)  = ^ ] x(iAT)  exp  (-jw  IaT) 


(3-16) 


The  time  domain  multiplication  of  the  original  signal  x(t)  by  the  impulse  train 
requires  that  the  Fourier  transform  of  the  resulting  sampled-data  sequence  be 
the  convolution  of  the  original  data  transform  with  the  transform  of  the  impulse 
train.  The  Fourier  transform  Is  readily  seen  to  be  a sequence  of  singularities 

at  frequencies  w = n||-»  1 = o,  ±1,±2 This  property  readily  follows  from 

the  fact  that  6y(t)  Is  periodic  with  periodicity  aT.  Also 


ao 

^ C<5j(t)J  = ^ ' exp  (-jwIaT) 


= <*>  for  w = 


(3-17) 


We  shall  show  that  the  sampled-data  Fourier  transform  is  related  to  the  original 
signal  transform!  by  the  equation 


oo 

**<“>  ■ ST  £ *(»  + if  ") 


(3-18) 
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where 


X(w)  = F[x(t)] 


Typical  plots  of  |X*(w)|  given  |X(w)|  are  plotted  in  figure  22.  The  extremely 


1X(»)| 
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Figure  22.  Spectral  Folding  of  Sampled  Signals. 

Important  result  is  that  frequency  components  in  X(w)  above  the  Nyquist  frequency 
w = tt/aT  are  folded-down  to  frequencies  within  the  band  |w|  ^ tt/aT.  Thus  these 
higher  frequency  components,  mimic,  or  alias  spectral  components  in  the  lower 
frequency  band.  A time  domain  interpretation  of  the  aliasing  phenomena  is  shown 
In  figure  23.  Observe  the  spectral  distortion  which  occurs  in  the  case  that  the 
sample  rate  Is  too  low.  Distortion  Is  minimized  with  sample  frequencies  at  least 
twice  the  signal  bandwidth.  Signals  which  are  bandlimited  to  the  Nyquist  fre- 
quency can  be  uniquely  represented  by  their  sample  sequence  x^ . Unfortunately, 
perfectly  bandlimited  signals  are  seldom  encountered  in  practice.  Many  practical 
signals  are  approximately  bandlimited,  however.  Those  signals  not  sufficiently 
bandlimited  for  the  desired  sample  rate  can  be  appropriately  filtered  by  anti- 
aliasing filters  to  minimize  their  spectral  content  beyond  the  Nyquist  frequency 
thereby  minimizing  aliasing  distortion.  Anti-aliasing  filtering  must  be  accom- 
plished prior  to  digitization. 
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Figure  23.  A Time  Domain  Example  of  Frequency  Aliasing. 

The  Important  relationship  of  equation  (3-18)  is  proven  In  the  following 
manner.  First  note  that  X*(w)  as  expressed  in  equation  (3-15)  is  periodic  in  w 
with  periodicity  wp  = 2tt/aT.  This  periodicity  is  a consequence  of  the  periodic- 
ity of  the  complex  exponential  function. 


exp  [-j  (|y  n+w)  ATl]  = exp(-j2irni)  exp  (-j  wATi) 

= exp  (-jwATi)  (3-19) 

since  exp(-j2iTni)  = 1 for  all  Integers  n and  i.  The  periodic  nature  of  X*(w)  per 
mits  it  to  be  represented  by  a sum  of  sinusoids  (Fourier  series)  of  the  form 


X*(w) 


oo 


c ( 1 ) exp  (-j'waTI) 
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where  the  Fourier  coefficients  c(i)  are  given  by 


_+w 


1 


P/2 


X*(w)  exp 


P/2 


[•»(?')] 


dw 


Comparing  equations  (3-15)  and  (3-20)  we  see  that 


c(1)  = xi  = x ( i AT ) 


and  equation  (3-21)  can  be  re-expressed  in  the  form 


x(iAT)  = 


u/AT 


X*(w)  exp  (+jwATi)  dw 


ir  / AT 


(3-21) 


(3-22) 


(3-23) 


Standard  Fourier  analysis  of  the  original  signal  expressed  as  in  inverse  trans- 
form is 


■fe/ 


x(t)  = y / X(w)  exp  (jwt)  dw 


(2m+l  )tt 
°°  /*  AT 

■tEf 


X(w)  exp  (+jwt)dw 


m=-oo  ) it 

AT 


AT 


/ 1 i 

2-r  X (w  + aT  m)  exp  l) 

ir  m=-« 


dw 


AT 


(3-24) 


;■ 
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I 


\ 


Restricting  attention  to  values  of  t = iAT  In  the  above  equation  and  making  the 
simpl if icaticn  in  the  exponential  that  this  allows  gives 


exp  (jwATi)  dw 


(3-25) 


Comparing  equations  (3-23)  and  (3-25)  and  the  uniqueness  of  the  Fourier  transform 
gives  the  identity 


x*w  = y Y, 

m=-“ 

as  was  to  be  proven. 

Up  until  this  point  we  have  assumed  that  an  infinite  duration  of  data  is 
available.  The  practical  considerations  of  a finite  duration  data  interval  dis- 
cussed In  section  III. 2 also  hold  for  sampled-data.  This  finite  data  restriction 
reduces  equation  (3-15)  to  the  form 

N-l 

X*T(w)  = y ' x.  exp (-jwaTI) 

i=0  (3-27) 


The  above  representation  of  the  Fourier  transform  of  the  signal  x(t)  is  simply 
computed  for  a given  w as  the  weighted  sum  of  exponentials.  We  shall  show  that 
X*^.(w)  need  be  evaluated  for  only  N equally  spaced  values  of  w in  order  to 
uniquely  specify  the  sampled  sequence  x...  Then  we  shall  show  how  these  N values 
may  be  efficiently  calculated. 

The  discrete  Fourier  transform  (DFT)  is  defined  simply  as  the  Fourier  trans- 
form of  a finite  length  sample  sequence  evaluated  at  the  discrete  frequencies 
w = n^j-,  n » o,  1,  2 , . . N • 1.  We  now  show  the  correspondence  between  the  DFT 
and  the  transform  X(w)  of  the  original  unsampled  signal. 
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Consider  the  windowed  time  function  x(t)  defined  on  the  interval  [0,T)  and 
elsewhere  zero.  Let  x^,  1 s o,  1,  . . . ,N-1  be  the  corresponding  sample  se- 
quence taken  at  times  t = iAT,  where  aT  = T/N.  Since  x(t)  is  zero  outside  the 
interval  [0,T),  we  may  assign  values  to  a function  x(t)  such  that  x(t)  = x(t)  for  , 

t in  [0,T)  and  in  such  a manner  that  allows  us  to  most  easily  reconstruct  x(t) 
from  x(t)  and  its  transform.  One  potentially  useful  assignment  is  to  define 
x(t)  as  the  periodic  repetition  of  x(t)  as  shown  in  figure  24. 


Figure  24.  Periodic  Continuation  of  a 
Finite  Duration  Signal. 


The  advantage  of  this  choice  is  that  because  of  its  periodicity,  it  can  be  rep- 

o 

resented  by  a Fourier  series  at  frequencies  wp  = -y  and  its  harmonics.  Further, 
if  x(t)  is  truly  bandlimited  at  wl  = then  on^  \/wp  = N/2  components  need  be 
calculated  to  uniquely  characterize  x(t).  We  shall  show  that  knowledge  of  N DFT 
points  is  sufficient  to  characterize  x^.  Knowledge  of  transform  values  at  other 
frequencies  is  not  required  and  in  fact  these  values  may  be  deduced  from  the  dis- 
crete frequency  transform  values  that  are  actually  evaluated. 

In  terms  of  a complex  Fourier  series,  x(t)  can  be  written 
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x(t)  = ^2  Xi  (k)  exp  f+j|2-  kt) 
k=-~  ' 7 


(3-29) 


or  when  t is  restricted  to  values  iAT 


oo 

c(iAt)  - XJ  (k)  exp  f^kATi) 


£ Xj  (k)  exp  (j^  ki) 
k»—  ' 


(3-30) 


since  T/aT  = N.  The  exponential  term  is  periodic  in  k with  periodicity  N 'so  that 
equation  (3-30)  may  be  written  as  a finite  sum 


*(UT>  ■ fr  £ 


Xp(k)  exp  kij 


(3-31) 


The  transform  Xp(k)  is  seen  to  be 


Xp(k)  - N ^ XJ  (k+nn) 


(3-32) 


The  transform  Xp(k)  is  most  readily  evaluated  from  equation  (3-31 > rather 
than  by  use  of  equation  (3-32).  Multiplying  both  sides  of  equation  (3-31)  by 
exp  (-j  jp-  it)  and  summing  over  the  index  1 gives 
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The  term  within  the  brackets  Is  recognized  as  a geometric  series  which  Is  sum- 
mable  In  closed  form.  It  equals  1_  whenever  (k  - i)  is  an  Integer  multiple  of  N 
and  is  zero  otherwise.  Hence  we  have  shown  the  transform  pair 


x(1aT)  <=>  Xp(k) 

The  transform  relationships  are  summarized  by  the  equations 


Xp(k)  = x(iAT)  exp  ^-j  1k^ 

i=o 

and 

N-l  . 

x ( i AT ) = J-  ^ Xp(k)  exp  ki  j 

k=o 


(3-35) 


(3-36) 


Equations  (3-35)  and  (3-36)  are  the  DFT  and  inverse  DFT,  respectively,  of  the 
sampled-data  sequence  x^  = x(iAT).  By  reconstructing  the  steps  leading  to  equa- 
tion (3-35),  the  relationship  between  Xp(k)  and  X(w),  the  Fourier  transform  of 
the  original  data  signal,  is  established.  These  steps  are  summarized  in  table  2. 
Observe  that  the  order  in  which  Steps  1 and  2 are  executed  may  be  reversed  with- 
out changing  the  result.  The  important  conclusion  to  be  drawn  from  table  2,  given 
the  assumptions  that 
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Table  2 

THE  RELATIONSHIP  OF  DISCRETE  FOURIER  TRANSFORMS  TO  THE 
ORIGINAL  SIGNAL  TRANSFORM 


Process  Input 


Output 


Transform  Relationshii 


ao 

Sampling  x(t)  x*(t)  = x(t)*6y(t)  X*(w)  = ^ * X (w+aT)1) 


2 Windowing  x*(t)  x|(t)  = g(t) -x*(t)  X|  = X*(w)*  G(w) 


3 Discrete  x. 
Fourier  1 
Transform 


x j = xf(iAT) 


Xp(k)  ■ X*(w)  l^l* 


(1)  the  sampling  rate  is  sufficiently  fast  that 


X*(w)  * ^ X(w)  , |w|  < ^ 


(2)  the  data  analysis  interval  is  sufficiently  long  that  the  frequency  dis- 
tortion induced  by  windowing  is  negligible  (even  though  a correction 
factor  may  be  required). 


is  that 


x (k)»  l„  = k 2, 

P aT^CF  W K T 


•w-k  ^»AT  Xp(k)  F 


(3-38) 
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where  Xp ( k ) Is  the  DFT  of  data  sequence.  It  is  defined  by  equation  (3-35).  This 
method  for  computing  X(w)  and  the  definitions  of  power  and  cross  spectral  densi- 


ties, transfer  function  estimates  and  the  coherence  function  establish  the  approx- 
imations 


PSD 


*xx  (r)  = T lXT(w)  I2 


*2JZ|X(k)|2T 

w=f  k N2  P 


CSD 


/K\  ^j(w)  V*(w) 

V ft) ' 2 T r- 


,CF 


w=k  p- 


2 , XD(k)YD('k)T 

N2  P P 


TRANSFER  FUNCTION  ESTIMATES 


(Method  1 ) 


H"  (k  r) = V ft)'  *xx  (f) 


■ Vk)'  Vk> 


(Method  2) 


I"  (k  r)  l * *1'/  ft)  »;)/4ft) 


|Yp(k)/  xp(k)| 


Coherence  Function 


'xy 


ft) 


1 *; 

<y_ 

ft) 

1 I2 

4xx  1 

ft) 

$ | 
yy  ’ 

ft) 

= 1 (see  note) 


(3-39) 


(3-40) 


(3-41) 


(3-42) 


(3-43) 


Note:  The  coherence  function  as  evaluated  here  Is  a highly  biased  estimate  of 
the  actual  true  value  of  the  coherence.  This  difficulty  Is  removed  by 
application  of  frequency  averaging  techniques  developed  in  section  III. 4. 
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Recall  that  the  above  equations  represent  estimated  quantities  and  the  above 
expressions  are  "raw"  estimates.  Averaging  techniques  for  reducing  the  variance 
of  the  estimates  are  developed  in  section  III. 4. 

The  use  of  equations  (3-39)  through  (3-43)  became  particularly  attractive 
for  digital  computer  evaluation  of  signal  properties  with  the  Introduction  of  the 
Fast  Fourier  Transform,  FFT,  by  Cooley  and  Tukey  (reference  9)  In  1965.  The  FFT 
Is  a computational  efficient  algorithm  for  evaluation  of  the  DFT,  equation  (3-35). 

The  FFT  algorithm  depends  upon  two  key  aspects.  One  aspect  is  that  the  com- 
putational time  required  to  evaluate  Xp(k),  k = 0,  1 , . . . N-l  Is  dominated  by 
the  time  required  to  compute  a complex  multiplication.  An  N point  sequence  re- 
quires N2  multiplications,  a computational  burden  which  Is  prohibitive  for  large 
N.  The  FFT  'trades'  complex  multiplications  for  an  increased  number  of  computa- 
tionally less  demanding  complex  additions  to  achieve  a computational  complexity 
comparable  to  N log2N  complex  multiplications.  The  computational  savings  is  the 
ratio 

Al2  _ N 

NTogp'  " iog2N  (3_44) 

The  computational  savings  ratio  is  tabulated  in  table  3 for  typically  used  values 
of  N. 

Table  3 

THE  COMPUTATIONAL  SAVINGS  RATIO  OF  THE  FAS  I FOURIER  TRANSFORM 


N 

N/log2N 

512 

57 

1024 

102 

4096 

341 

16384 

1170 

For  example,  the  numerical  evaluation  of  a 512  point  DFT  by  straight  forward 
evaluation  of  equation  (3-35)  would  require  57  times  as  much  computer  time  as 
would  the  FFT  to  evaluate  the  same  512  point  sequence.  A second  aspect  of  the 
DFT  formula  which  makes  the  FFT  possible  Is  the  periodicity  of  the  complex 
exponential . 
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exp[  j jjp  (1+in)]  = exp  (j^-  i)  (3_45) 

E 

These  aspects  of  computing  DFT's  lead  to  computational  advantage  as  follows: 

First  we  observe  that  the  mathematical  expression  E = ab  + ac  + ad  can  be  evalu- 
ated explicitly  as  written  to  require  three  multiplications  and  two  additions. 

Alternatively,  E = a(b+c+d)  evaluated  explicitly  requires  two  additions  but  now 
j requires  only  a single  multiplication  to  obtain  an  equivalent  result.  This  con- 

cept Is  similarly  applied  toward  evaluation  of  the  DFT. 

We  separate  a N point  sequence,  assuming  N is  an  even  number,  into  two  N/2 
point  sequence  x2^  and  x2jl+j . One  sequence  has  been  formed  from  the  points  having 


even  indices  and  the  other  from  the 

points  having  odd  indices.  Applying  equation 

(3-35) 

N-l 

i \ 

DFT  (Xl)  - Xp(k)  - 

M 

X 

exp(_J’  r ik) 

1*0 

N/2-1 

■ E 1 

x2i  exp 

(-j  ^L21k) + x2i+i  exp  ("j  r 21+lk) 

I 

1 

1*0 

N/2-1 

■ E 

x21+x21+l  exp  (-j  IT ) exp  (-j  Wl  ik) 

1*0 

■ DfT  (x2l' 

| + exp  | 

■irl)1FT.(*2w) 

(3-46) 

Thus  the  original  N point  transform  has  been  changed  to  an  equivalent  problem  of 
evaluating  two  N/2  point  sequences  and  linearily  combining  the  results.  The 
computational  saving  accrues  because  the  computation  burden  is  proportional  to 
the  sequence  length  squared.  The  N2  multiplications  required  for  direct  evalua- 
tion of  the  DFT  is  reduced  to  2(N/2)2  + N»N2/2  multiplications. 

A secondary  advantage  of  the  FFT  technique  is  that  one  evaluation  of  the  com- 
plex exponential  can  be  used  for  the  evaluation  of  both  the  two  smaller  sequences. 
Thus,  fewer  evaluations  of  the  complex  exponential  are  required. 
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The  process  of  dividing  evaluation  of  an  DFT  Into  the  evaluation  of  two 
OFT's  having  half  as  many  points,  as  done  in  equation  (3-46),  can  be  repeated  as 
long  as  the  number  of  points  in  the  DFT  actually  evaluated  is  even.  Cooley  and 
Tukey  (reference  9)  developed  the  algorithm  for  special  case  that  N = 2m,  m an 
integer.  In  this  case,  N multiplications  are  performed  for  each  of  the  m steps 
in  which  DFT's  are  expressed  as  the  linear  combinations  of  two  shorter  DFT's. 

The  total  number  of  multiplications  required  by  this  method  is 

Nm  = nlog2N  (3_47) 

compared  with  N2  multiplications  required  for  direct  evaluation  of  the  DFT  given 
by  equation  (3-35). 

Other  methods  applicable  to  cases  in  which  N is  not  a power  of  2 have  also 
been  developed  (references  10,  11,  and  12).  These  methods  involve  the  decomposi- 
tion of  N into  its  component  prime  factors.  Each  of  these  methods  is  equally 
applicable  to  evaluation  of  the  inverse  DFT,  equation  (3-36). 

3 CONFIDENCE  INTERVALS 

In  section  II. 5 we  showed  that  practical  considerations  of  random  data  anal- 
ysis prevented  us  from  precisely  calculating  the  statistics  of  an  observed  pro- 
cess. Rather,  sample  statistics  can  be  calculated  which  approximate  the  desired 
process  statistics.  In  this  section  we  shall  derive  a measure  of  the  error  or 
quality  of  our  estimates. 

Generally,  the  measurement  error  by  which  a sample  statistic  differs  from 
the  true  quantity  is  a random  variable  since  it  is  the  function  of  random  vari- 
ables. Consequently,  mean  values,  moments  and  probability  density  functions 
characterize  the  error  random  variable.  This  aspect  of  sample  statistics  is 
illustrated  with  an  example. 

Given  independent  random  variables  x^  with  means  u and  variances  a2,  we  make 
N observations.  From  these  observations  we  wish  to  make  an  estimate  p of  the 
process  mean.  A logical  choice  is  to  simply  average  the  observations  to  form  the 
estimate. 


(3-48) 
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The  uncertainty  of  our  error  depends  upon  the  probability  of  the  underlying 
process.  Observe  that  for  a very  large  N,  the  mean  estimate  becomes  a very  good 
one  since  the  error  variance  decreases  as  N"1. 

The  estimate  error  Is  a random  variable.  The  question  of  exactly  how  much 
in  error  our  particular  estimate  is  cannot  be  answered.  In  many  cases,  however, 
we  can  compute  the  probability  that  the  error  exceeds  a specified  bound.  Both 
the  probability  and  bound  must  be  specified  to  characterize  the  'goodness'  of 
our  estimate. 

These  concepts  are  applied  to  the  mean  estimation  problem  developed  in 
equations  (3-48)  through  (3-50).  We  assume  that  the  random  variables  x^  have 
normal  distributions.  Then  the  estimated  mean  is  also  norma’ly  distributed  and 
has  a PDF  completely  characterized  by  its  mean  and  covariance  given  in  equations 
(3-49)  and  (3-50). 
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. 1 (p-u)2 
2 si 

N 


The  estimate  error  e = u-v  has  probability  density  function 


(3-51) 


VF 


1 e2“ 
exp  -2“ 
oz 

L rJ 


(3-52) 


By  the  definition  of  probabilities  presented  in  Section  II,  the  probability 
that  the  absolute  error  exceeds  a specified  bound  ed  is  precisely  calculable  from 
the  error  probability  density  function  equation  (3-52).  Explicitly, 


P ■/..  p(e)de  + J"  p(e)  de 


2 / p(e)  de 


(3-53) 


Table  4 gives  evaluations  of  this  Inteqral  for  selected  values  of  normalized 

Table  4 

CONFIDENCE  PROBABILITY  TABLE, a = aN’l/2 


ed/CTe 


P({|e|  > ed}) 


P({|e|  < e.}) 
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by  CTg^crN"1/ 2 . Given  = 2ag,  for  example,  we  say  that  there  is  only  a 4.6% 
probability  that  our  estimate  error  magnitude  exceeds  2 afi.  Equivalently,  we  say 
that  we  have  95.4%  confidence  that  our  estimate  error  magnitude  is  less  than  2ag. 

An  Important  aspect  of  the  preceding  analysis  if  that  two  quantities  must  be 
enumerated  to  Indicate  the  quality  of  the  estimate.  The  required  quantities  are 
(1)  the  error  bound  and  (2)  the  probability  that  the  specified  bound  Is  (or,  is 
not)  exceeded.  Generally  one  of  several  alternative  approaches  is  used.  Calcu- 

s 

latlon  of  the  confidence  from  an  a priori  error  bound  is  one  approach  commonly 
used.  Otherwise,  a particular  confidence  is  specified  and  the  corresponding 
error  bound  is  calculated.  An  approach  we  shall  use  entails  specifying  both  the 
error  bound  and  the  confidence.  Once  these  quantities  are  specified,  equations 
(3-52)  and  (3-53)  are  used  to  calculate  a minimum  number  of  observations  N that 
must  be  made  to  obtain  the  desired  quality  estimate. 

An  identical  analysis  procedure  is  applicable  to  the  calculation  of  confi- 
dence bounds  for  power  spectral  densities  and  cross  spectral  densities.  Results 
presented  in  references  2,  8,  and  12  are  summarized. 

The  underlying  probability  of  measured  signals  of  dynamic  systems  tends  to 
be  a normal  distribution  as  demonstrated  by  the  central -limit  theorem.  It  fol- 
lows that  the  real  and  imaginary  parts  of  a signal  DFT  at  each  component  fre- 
quency are  independent,  normally  distributed  random  variables  with  zero  means 
and  equal  variances  (reference  2,  p.  189).  The  PSD  is  calculated  from  the  real 
and  imaginary  DFT  components  by  equation  (3-39)  or 

iVk)'JT 

(3-54) 

As  the  sum  of  squares  of  normal  random  variables,  $xx  (y-)  is  a chi-square  random 
variable  with  two  degrees  of  freedom  (reference  4,  p.  250).  The  ratio  of  the 
standard  deviation  to  the  mean  of  a. chi-square  distribution  having  n degrees  of 
freedom  is  a constant  given  by 
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Chi-square  probability  density  functions  for  higher  degree  of  freedom  are 
plotted  in  figure  25. 


Figure  25.  Chi-Square  Probability  Density  Functions  for 
Several  Degrees  of  Freedom. 

Clearly,  our  PSD  estimate  Is  a chi-square  random  variable  with  two  degrees 
of  freedom.  As  a consequence  of  equation  (3-55)  with  n equal  to  2,  the  PSD  e<- 
mate  given  by  equation  (3-54)  has  a variance  equal  in  magnitude  to  the  est*-  • 
Itself.  This  estimation  error  Is  Intolerably  large  for  most  appl icatic''' 

♦Roughly  speaking,  the  number  of  degrees  of  freedom  is  the  re- 
normal  random  variables  summed  to  form  the  chi-square  variat'- 

That  Is  x2  = xx2  + X,2  + . . . x 2 has  n degrees  of  free:  ^ 

(reference  4,  p.  250). 
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The  estimation  errors  of  the  raw  PSDs  computed  by  equation  (3-54)  are  re- 
duced by  averaging  techniques.  Averaging,  In  effect.  Increases  the  number  of 
degrees  of  freedom  of  the  chi-square  random  variable  or  PSD.  By  equation  (3-55) 
this  averaging  process  reduces  the  estimated  variance  In  proportion  to  the  esti- 
mate mean,  a constant.  The  variance  reduction  Is  proportional  to  *“1/2  where  i 
Is  the  number  of  Independent  raw  PSDs  averaged  to  obtain  the  estimate. 

Two  distinct  averaging  techniques  are  commonly  employed.  One  technique  Is 
to  average  several  PSD  estimates  each  obtained  from  different  time  Intervals  of 
the  data  record.  The  second  technique  Is  to  average  several  adjacent  frequency 
components  of  the  original  raw  PSD. 

Averaging  over  several  PSD  estimates  (smoothing  over  an  ensemble  of  mea- 
surements) Increases  the  number  of  degrees  of  freedom  of  the  chi-square  dis- 
tributed estimate  as  Is  seen  from  the  following  calculations.  Denoting  the 

A 

averaged  estimate  by  we  have 


[ 

K 

' ' 

r i 


xx 


N 1=1 


{v  [Vk>] 


Im2  [XF 


,(k)]} 


(3-56) 


Then  If  i raw  PSD  estimates  are  averaged,  the  smoothed  PSD  estimate  has  n * 2* 
degrees  of  freedom.  Application  of  equation  (3-55)  to  this  average  shows  that 
the  standard  deviation  of  the  smoothed  estimate  Is  i~1/2  times  less  than  that  of 
the  original  estimate. 

Similar  results  are  obtained  for  smoothed  estimates  derived  as  averages  of 
adjacent  frequency  components  of  a single  raw  PSD.  This  averaging  technique  Is 
particularly  advantageous  whenever  a *xx  (or  log  *xx)  versus  log-frequency  pres- 
entation of  the  PSD  Is  employed  as  Is  commonly  the  case  (reference  8,  pp.  90-93). 
The  spacing  of  PSD  estimates  uniformly  In  frequency  as  occurs  with  DFT  computa- 
tions gives  a far  greater  resolution  that  can  be  leglblly  displayed  on  a loga- 
rlthmetlc  frequency  graph.  These  'extra'  PSD  estimates  can  be  averaged  to  reduce 
the  uncertainty  of  the  estimate.  For  frequency  averaging. 
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(3-57) 


An  apparent  disadvantage  of  the  frequency  smoothing  technique  In  comparison 
with  the  ensemble  smoothing  technique  of  equation  (3-56)  Is  the  loss  of  obtain- 
able resolution,  especially  at  low  frequency  ranges.  However,  the  ensemble 
smoothing  technique  requires  additional  data  samples.  For  an  Identical  statis- 
tical confidence  of  PSD  estimates  given  the  same  number  of  signal  data  points, 
each  method  yields  the  same  frequency  resolution.  The  ensemble  smoothing  tech- 
nique Is  computationally  less  expensive,  however.  If  N points  are  divided  Into 
i Intervals  for  which  transforms  and  PSD's  are  computed  and  averaged,  then  the 
computational  savings  over  frequency  smoothing  a N point  sequence  to  obtain  the 
same  confidence  Is  roughly  a factor  log2t. 

That  the  ensemble  averaging  and  frequency  smoothing  yield  the  same  confi- 
dence and  frequency  resolution  given  the  same  number  of  raw  data  points  Is 
Illustrated  by  the  following  example. 

EXAMPLE  3-1 

Given  N = 4 x 210  data  points,  compute  the  frequency  resolution  of  smoothed 
PSD  estimates  having  a confidence  equivalent  to  8 degrees  of  freedom.  A 1 KHz 
sampling  frequency  Is  used. 

Method  1,  Ensemble  smoothing: 

Four  separate  DFTs  must  be  computed  to  yield  the  required  degree  of  freedom.  One 
fourth  the  total  available  samples,  or  210  points,  is  available  for  each  DFT. 

The  frequency  resolution  obtained  Is 


Af  ■ i " — “ — - 1.0  Hz 
T 21°/fs 


Method  2,  Frequency  smoothing: 

Four  adjacent  PSD  frequency  components  must  be  averaged  to  yield  an  estimate 
with  the  required  degrees  of  freedom.  In  this  case,  all  N data  points  are 
available  for  the  raw  DFT  computation.  The  unsmooth  frequency  resolution  Is 
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This  frequency  resolution  Is  4 times  smaller  than  available  In  Method  1 since  the 
data  Interval  used  for  the  DFT  calculation  Is  4 times  longer.  But  we  must  aver- 
age 4 adjacent  frequency  components  to  obtain  the  desired  confidence. 

if  * 4 • Af 1 = 1 . 0 Hz 

Identical  frequency  resolution  Is  obtained  In  either  case.*-* 

Confidence  bounds  are  computed  from  a normalized  chi-square  distribution. 
Normalizing  this  random  variable  by  Its  mean  we  obtain  multiplying  factors  which 
determine  the  upper  and  lower  confidence  bounds  of  PSDs  as  a multiple  of  the 
estimate  Itself.  Small  confidence  bounds  require  that  the  multiplying  factor  be 
nearly  equal  to  unity. 

Upper  and  lower  confidence  bound  multiplying  factors  tabulated  In  table  5 


Table  5 

CONFIDENCE  INTERVAL  UPPER  AND  LOWER  BOUNDARY  MULTIPLIERS, 
M,  FOR  AN  n DEGREE  OF  FREEDOM  CHI-SQUARE  DISTRIBUTION 


100a: 

80  Percent 

90  Percent 

95 

Percent 

21s  n 

M: 

Mu 

Mu 

Mu 

2 

.n 

2.31 

.05 

3.00 

.03 

3.69 

10 

.49 

1.60 

.39 

1.83 

.33 

2.05 

20 

.62 

1.42 

.54 

1.57 

.48 

1.71 

60 

.77 

1.24 

.72 

1.32 

.68 

1.39 

120 

.83 

1.17  ' 

.80 

1.22 

.76 

1.27 

are  derived  from  the  normalized  chi-square  distribution.  Both  upper  and  lower 
factors  must  be  specified  since  the  chi-square  density  function  Is  asymetrlc 
about  Its  mean.  Calculation  of  upper  and  lower  multiplying  factors  for  cases 
not  presented  In  table  5 are  calculable  by  methods  Illustrated  In  Example  3-2. 

EXAMPLE  3-2 


Given  the  percent  confidence  100  a that  we  desire  a PSD  estimate  to  have  and  the 
number  of  degrees  of  freedom  n equal  to  twice  the  number  of  raw  DFT  components 
averaged  to  calculate  the  PSD  estimate,  calculate  $he  upper  and  lower  confidence 
bound  multipliers  for  the  estimate.  The  estimate  has  a chi-square  distribu- 
tion with  mean  * and  standard  deviation  [see  equation  (3-55)] 
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-V§ 


where  n Is  the  number  of  degrees  of  freedom  and  Is  equal  to  7JL.  Thus  we  need  to 
calculate  the  upper  and  lower  multiplier  Mu  and  such  that 


(fxx  ” H„  *xx})  + P({*xx  * Ht  *xx}) 


To  ensure  an  unbiased  confidence  Interval  we  further  require  that 


P 4xx  » H 


'»  fxx})  • P({ 


$ < M.  ♦ : 

XX  l XXI 


a new  random  variable  y * n #xx^®xx  defined  to  facilitate  later  use  of  mathe- 
matical tables.  Certainly,  y nas  chi-square  distribution  with  n degrees  of  free- 
dom and  mean  n since  It  Is  simply  a scalar  multiple  of  $xx.  Now  the  multiples 
can  be  formulated  as 


P ({‘xx  * Mu  \xl)  ■ p({*  - "M) 


' f p(y)  dy  ■ \ (l-o) 


where  p(y)  Is  the  n degree  of  freedom  chi-square  probability  density  function. 

Tabulation  of  nM  given  n and  the  Integral  value  have  been  published  (ref- 
erence 13).  u 

For  the  purpose  of  Illustration,  let  o * .9  corresponding  to  90*  confidence 
and  suppose  n ■ 10.  Then  from  the  tables  we  find  that: 


p(y)  dy  - .05 


for  nM  ■ 18.3  .corresponding  to  Mu  - 1.83.  To  compute  the  lower  multiple  M^ 
using  'ihe  tables,  the  problem  must  be  recast  Into  the  form 
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p(y)  dy  * [ p(y)  dy  ■ i - \ (l-o) 


p(y)  dy  - .95 


From  the  tables,  nMa  a 3.9  or  M$,  a 0.39.  These  values  agree  with  those  tabu- 
lated In  table  5.  An  Identical  calculation  technique  Is  applied  when  evaluating 
other  cases.*-* 

Often  PSD  estimates  are  presented  In  a log  magnitude  versus  log-frequency 
graphical  format.  For  this  presentation,  a confidence  Interval  'spread'  ex- 
pressed In  decibels  Is  more  convenient  to  use  than  the  multipliers  and  Mu- 
Taking  logarithms 


S i 10  log10  ^ 


10  1 og x o CMU ) - 10  log10(M4) 


(3-58) 


The  spread  S expressed  In  dB  Is  tabulated  In  table  6 as  a function  of  confidence 


Table  6 

CONFIDENCE  INTERVAL  SPREAD  (dB)  FOR  AN  n DEGREE 
OF  FREEDOM  CHI-SQUARE  DISTRIBUTION 


100a: 


80  Percent 
13.4 


90  Percent 
17.6 

6.7 

4.6 

2.6 

1.8 


95  Percent 
21.7 
11.9 
5.5 
3.1 
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and  the  chi-square  degrees  of  freedom  of  the  estimate.  Blackman  (reference  14, 
p.  23)  presents  a fairly  accurate  formula  for  calculating  the  dB  confidence 
spread  associated  with  chi-square  distributions.  The  formula,  reasonably  accu- 
rate for  n > 10  Is 


(3-59) 


where  K Is  determined  by  the  desired  confidence  according  to  table  7.  For  a 


Table  7 


CONFIDENCE  SPREAD  CONSTANTS 
(See  equation  (3-59)) 


K(dB) 

16 


V 


given  confidence,  S Is  readily  calculated  as  a function  of  n.  Conversely,  n can 
be  expressed  as  a function  of  S;  this  technique  Is  commonly  employed  to  deter- 
mine the  number  of  terms  that  must  be  averaged  to  achieve  a specified  confidence 
and  error  tolerance  bound. 

Similar  concepts  may  be  applied  toward  determining  confidence  bounds  for 
CSD's,  coherence  functions,  and  transfer  function  estimates.  The  procedures  are 
significantly  more  difficult  In  these  cases.  The  references  present  some  of  the 
details  (reference  2,  pp.  193-203).  These  details  are  not  presented  In  this 
guide.  Analyses  for  these  cases  show,  as  In  the  case  of  PSD  estimates,  that  the 
quality  of  the  estimate  Improves  as  more  "raw"  estimates  are  averaged  to  form  the 
best  estimate.  Averages  estimates  for  CSD's  should  be  formulated  In  the  same 
manner  as  averages  for  PSD's.  Averages  for  coherence  functions  and  transfer 
function  estimates,  however,  should  be  computed  from  averaged  PSD's  and  CSD's 
rather  than  by  averaging  raw  coherence  and  transfer  function  estimates.  This  Is 
readily  seen  to  be  valid  In  the  case  of  coherence  function  estimates  where  the 
later  procedure  would  always  yield  estimates  Identically  equal  to  unity  (see 
equation  (3-43)),  which  Is  obviously  an  Incorrect  result.  Furthermore,  recall 
the  technical  development  of  section  II. 9 which  demonstrated  that  transfer 
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functions  estimates  are  only  valid  over  the  frequency  ranges  In  which  the  co- 
herence function  tends  toward  unity.  Otherwise,  the  "output"  Is  not  strongly  In- 
fluenced by  the  so-called  "Input"  signal  and  a transfer  function  concept  Is  In- 
appropriate. 

5.  SPECIAL  TOPICS 


Several  additional  topics  of  Importance  In  practical  data  analysis  are  dis- 
cussed briefly  In  this  section.  The  topics  discussed  are: 

• Anti-aliasing  filters 

• Special  considerations  for  CSD  estimates 

• Appending  zeroes  to  signal  sequences 

• Bias  and  trend  removal 

• Sinusoidal  signal  components 

• Transfer  function  estimates  for  deterministic  inputs. 

a.  Anti-aliasing  Filters 

Anti-aliasing  filters  are  required  to  prevent  (or  at  least  minimize)  the 
aliasing  or  mimicking  of  low  frequency  signal  components  In  a digitized  signal  by 
high  frequency  components  In  the  original  analog  signal.  The  filtering  must  be 
accomplished  prior  to  signal  digitization.  The  filter  attenuation  beyond  the 
Nyqulst  frequency  (one-half  the  sampling  frequency)  must  be  sufficiently  great 
that  the  amplitude  of  the  attenuated  high  frequency  component  Is  small  in  compari 
son  with  the  magnitude  of  the  corresponding  In-band  frequency  component  so  that 
aliasing  distortion  of  a signal  spectrum  Is  negligible.  This  requirement  dic- 
tates that  multipole  filters  be  used  to  give  adequate  out-of-band  attenuation. 
Also,  we  require  that  these  anti-aliasing  filters  have  a minimum  distortion  ef- 
fect on  the  signal  spectrum  for  a reasonably  large  portion  of  the  useable  fre- 
quency band.  For  these  reasons,  "sharp  cutoff"  filters  are  used  having  break 
frequencies  equal  to  approximately  40%  of  the  sampling  frequency  f$.  Signal  fre- 
quency components  below  .4f$  are  passed  with  negligible  attenuation  while  compo- 
nents above  .4f$  are  attenuated  several  orders  of  magnitude.  Signal  phase,  on 
the  other  hand.  Is  noticeably  affected  throughout  the  filter  passband.  This  sig- 
nal phase  distortion  Is  inconsequential  for  PSD  estimates  but  becomes  critical  In 
CSD  and  transfer  function  estimates.  In  any  event,  the  data  analyst  must  under- 
stand the  purpose  of  the  anti-aliasing  filter  and  the  manner  in  which  It  affects 
the  Interpretation  of  data,  especially  If  non-standard  filters  are  chosen. 
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b.  Special  Considerations  for  CSD  Estimates 

Unlike  PSD  estimates,  CSD  estimates  require  calculations  Involving  two  dis- 
tinct signals.  The  CSD  estimate  Is  sensitive  to  any  process  which  affects  the 
relative  phase  of  the  two  signals.  This  phase  distortion  addltlvely  corrupts 
the  CSD  phase  estimate.  Special  attention,  then,  must  be  paid  to  sources  of 
relative  phase  distortion  whenever  the  CSD  phase  Information  Is  Important.  One 
source  of  relative  phase  distortion  Is  difference  In  signal  phase  due  to  use  of 
different  anti-aliasing  filters  having  differing  phase  properties.  Anti- 
aliasing filters  with  Identical  phase  characteristics  must  be  used  If  the  CSD 
phase  estimate  Is  not  to  be  distorted.  Alternatively,  the  relative  phase  dis- 
tortion between  the  filters  can  be  measured  or  calculated  and  that  result  used 
to  correct  the  CSD  phase  computation.  Another  primary  source  of  relative  phase 
distortion  between  two  signals  Is  time  delay  between  the  signals.  A pure  time 
delay  has  associated  with  It  a phase  delay  that  Is  linear  In  frequency.  The 
phase  difference  Is  related  to  the  time  delay  Td  by  the  equation: 


♦ =~wTd  = ~2irf  Td 


Generally,  this  phase  delay  Is  Intolerable  at  high  frequencies  for  Td  > (lOf  )-1. 
By  equation  (2-95),  this  phase  component  Is  Induced  Into  the  CSD  by  the  time  de- 
lay transfer  function.  Sources  of  an  equivalent  time  delay  between  two  signals 
Include  Inaccurate  time  bases  for  either  or  both  signals,  recording  delays  caused 
by  recorder  magnetic  head  skew,  electronic  recording  and  reproduce  delays,  and 
asynchronous  sampling  while  digitizing.  Whenever  accurate  phase  Information  is 
required,  these  time  delays  must  be  minimized,  and  the  residual  delay  must  be 
measured  so  that  It  may  be  removed  from  the  estimate  a posteriori. 

Another  consequence  of  signal  time  shift  Is  a reduction  of  apparent  coher- 
ence In  proportion  to  the  ratio  of  time  shift  to  the  data  Interval  time  duration. 
As  shown  In  reference  15: 


Y2  - y2  [1-t/TP 


where  t Is  the  time  shift  and  T the  data  Interval  duration. 
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c.  Appending  Zeros  to  Signal  Sequences 

Several  authors  (references  2 and  12)  recommend  appending  zeros  to  a data 
sequence  to  Increase  the  total  number  of  signal  points  to  an  Inteqer  of  two. 

Such  a total  number  of  points  Is  particularly  suitable  for  the  most  common  FFT 
algorithms  which  can  operate  only  upon  data  sequences  having  2m  points.  Adding 
zeros  effectively  alters  the  Intended  window  function  and  its  desired  spectral 
properties.  These  disadvantages  outweigh  any  advantage  gained  by  using  addi- 
tional data  points.  Alternatively,  modern  (reference  11,  p.  307)  FFT  algorithms 
have  been  devised  which  effectively  transform  N point  sequences  for  any  N that 
Is  highly  composite  (i.e.,  N Is  the  product  of  many  factors).  These  methods 
should  be  employed  In  cases  where  maximum  frequency  resolution  Is  required  and 
hence  the  longest  data  Interval  possible  must  be  transformed.  Also,  appending 
zeros  to  a data  sequence  will  change  the  RMS  level  calculated  for  the  signal 
Interval . 

d.  Trend  Removal 

For  many  applications,  the  signal's  mean  value  and  the  slope  of  a straight 
line  approximation  are  Important  signal  characterizations.  The  mean  value 
Identifies  dc  signal  bias  while  the  line  slope  Identifies  trends  or  a time  cor- 
relation of  the  data.  These  values  are  calculated  by  the  least-squares  tech- 
niques presented  In  section  II. 5.  Once  these  values  are  known,  no  additional 
Information  Is  gained  by  calculating  the  PSD  of  the  data  trend  line  component, 
as  the  corresponding  PSD  Is  uniquely  parameterized  by  the  dc  and  the  trend  slope. 
Trends  must  be  removed  If  measurement  of  the  low  frequency  random  spectra  Is  Im- 
portant to  the  analysis;  otherwise,  the  low  frequency  spectra  contributed  by  the 
trend  can  dominate  the  spectra  contributed  by  the  random  signal  component.  Win- 
dow functions  can  also  be  used  to  reduce  the  leakage  distortion  of  linear  trends 
(reference  8,  p.  81).  More  generally,  any  deterministic  signal  component  should 
be  removed  from  the  signal  prior  to  further  processing  to  characterize  signal 
stochastic  properties.  Also,  trend  removal  Is  Important  In  preventing  numerical 
Inaccuracies  that  may  arise  when  computing  transforms  on  small  word  length  com- 
puters . 

Procedures  for  removing  deterministic  components  from  measured  signals  must 
be  exercised  judiciously  so  that  only  actual  trends  are  removed.  Bendat  (refer- 
ence 2,  p.  291)  recommends  that  the  trends  be  removed  only  If  the  trend  Is  "phys- 
ically expected  or  clearly  apparent  In  the  data." 
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e.  Sinusoidal  Signal  Components 

Often,  spectral  techniques  Identify  resonances  or  modes  associated  with 
the  process  being  observed.  Lightly  damped  resonances  have  corresponding  PSD's 
exhibiting  large  spikes  at  one  or  more  discrete  frequencies.  These  spikes  are 
associated  with  sinusoidal  (or  periodic)  signal  components.  The  resonant  fre- 
quency and  the  energy  associated  with  the  sinusoid  are  typically  the  quantities 
of  most  Interest  to  the  data  analyst.  The  resonant  frequency  Is  determined  by 
the  frequency  at  which  the  PSD  has  a local  maximum.  The  uncertainty  of  the  fre- 
quency estimate  Is  the  resolution  bandwidth  of  the  PSD,  &f  equals  J./T  for  gen- 
eralized frequency  smoothed  PSD  estimates.  The  energy  associated  with  the 
resonance  Is  estimated  by  ascribing  all  the  energy  of  the  PSD  spike  to  the 
resonance.  This  energy  Is  obtained  by  the  techniques  developed  In  section  1 1. 9. 


E(fs) 


s + Af/2 
- af/2 


*xx(f)df 


(3-60) 


where  f Is  the  resonant  frequency  and  af  is  approximately  5 times  the  transform 
fundamental  frequency  resolution.  This  bandwidth  Is  required  to  Include  the 
spectral  leakage  side-lobe  components  as  shown  In  figure  20b. 

The  amplitude  of  a sinusoid  associated  with  the  spectral  resonance  spike  is 
estimated  from  the  spike  energy  with  the  assumption  that  all  the  spike  energy  is 
associated  with  the  energy  of  a sinusoid.  A sinusoid  with  ^ero-to- peak  amplitude 
A has  energy  ^A2.  It  follows  that 

A « y 2 • E (f$)  (3-6! ) 

Note  that  the  accuracy  of  the  amplitude  estimate  Is  proportional  to  the  square 
root  of  the  energy  estimate  error.  Thus  a 10%  energy  error  gives  a 5%  amplitude 
estimate  error.  Typically  the  energy  estimate  Is  high  since  It  Includes  energy 
from  the  random  signal  component  as  well  as  the  sinusoidal  component.  Good  ampli- 
tude estimates  are  achieved  provided  that  the  PSD  spike  is  an  order  of  magnitude, 
or  so,  higher  than  the  surrounding  PSD  floor,  as  shown  In  figure  26  for  a typical 
spike. 
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Figure  26.  A Sinusoidal  Signal  PSD  Spike. 


Other  methods  based  upon  the  absolute  PSD  spike  amplitude  are  available  for 
sinusoid  amplitude  estimation.  These  techniques  are  generally  less  accurate  than 
the  energy  method  because  the  spike  amplitude  depends  upon  the  window  function 
and  the  type  of  averaging  employed.  Additionally,  the  spike  amplitude  may  vary 
as  much  as  50%  depending  upon  the  alignment  of  the  true  sinusoid  frequency  with 
the  DFT  frequencies. 

f.  Transfer  Function  Estimates  for  Deterministic  Inputs 

One  has  a natural  tendency  toward  using  the  transform  technique  developed  in 
this  guide  for  determining  system  transfer  functions  given  samples  of  determinis- 
tic system  inputs  and  outputs.  In  view  of  the  stochastic  process  results,  we  are 
further  tempted  to  simply  ratio  output  to  Input  DFT's  to  obtain  system  transfer 
function  estimates  directly.  These  concepts  can  be  applied  provided  that  the 
analyst  bears  In  mind  the  limitations  of  the  technique  and  the  possible  gross 
estimate  distortions  that  may  result.  These  distortions  arise  from  the  aliasing 
and  leakage  distortion  associated  with  DFT's.  Transforms  computed  by  this  method 
correspond  to  systems  having  periodic  steady  state  waveforms  identical  to  the  ob- 
served transient  waveform.  Problems  associated  with  applying  transform  techniques 
to  deterministic  transfer  function  estimation  are  illustrated  in  the  following  ex- 
ample. 
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EXAMPLE  3-3 


Consider  the  system  in  figure  27a  excited  by  a step  input  initiated  at  time  T 


a.  System 


V JLI  T 

N 

b.  Transient  Input  and  Output 


c.  Periodic  Steady-State  Response 


Figure  27.  System  Identification  From  Deterministic  Input  Output  Data, 

The  process  is  observed  from  time  zero  until  T with  data  samples  taken  T/N  sec 
onds  apart.  The  output  is 


y(t)  ■ ^ e-3^  t^u_i(t-Ts 


fl  - e'^-Vl  u_,  (t-TJ 


t 
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where  the  function  u_:  (t-T  ) is  defined 


x(t)  * u_i  (t-T$) 


. f 1 ’t“T 

( 0 , otherwise 


the  corresponding  DFT's  are  computed  to  be 


/ v ^ / 2 it  \ exp  ("j  r* *k)  ' 1 

Vk)  3 ZL#  exp  (-J  1k)  7 2tt  T 

p Vi  N 1 -exP  (-jfk) 


and 


Yn(k) 


[l  - exp  (-a1?T/N)]exp  (-j  ^ Ik) 
W 


exp  | 

(-jjjp  *k)  - 1 exp 

>T  0 - w) 

j - exp 

[“  fr 

(j  2*k)' 

1 - exp  (-j  k) 

1 - exp 

[- 

1 (j  2 irk  + 

aT)‘ 

A transfer  function  estimate  formed  by  ratiolnq  Y (k)  and  X (k)  Is  compared  with 

(O^k \ -I  r P 

1 + j^j-)  In  figure  28.  Note  the  peculari- 


Itles  that  occur  In  even  this  most  elementary  example. 


It  follows  from  the  Fourier  series  interpretation  of  the  discrete  Fourier 
transform  that  a perfect  transfer  function  estimate  is  obtained  only  In  the  case 
that  the  transforms  of  steady  state  Input  and  output  waveforms  as  shown  In  figure 
27c  are  used  in  the  analysis,  and  that  the  input  signal  Is  bandl  Imlted.*-* 

As  Illustrated  In  the  example,  particular  care  must  be  exercised  when  ex- 
tracting function  estimates  from  deterministic  input  output  data.  Similar  care 
is  not  generally  required  with  ergodlc  random  data  for  which  each  analysis  inter- 
val exhibits  the  random  character  of  the  entire  data  interval. 
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SECTION  IV 

A NUMERICAL  DATA  ANALYSIS  EXAMPLE 


1 INTRODUCTION 

An  appreciation  of  the  benefits  and  utility  as  well  as  the  limitations  of 
random  data  analysis  techniques  can  only  be  gained  through  practical  experience 
In  the  use  of  these  techniques.  In  this  section,  a simplified  dynamic  system 
with  random  Inputs  Is  formulated  to  demonstrate  the  analysis  procedures  and  the 
proper  Interpretation  of  results.  Theoretical  analysis  (section  II  technique)  Is 
applied  to  the  system  model  to  predict  the  results  we  would  expect  to  obtain  by 
observing  system  Inputs  and  responses,  and  analyzing  these  data  with  practical 
data  analysis  procedures.  A numerical  model  of  the  system  Is  developed  and  simu- 
lated with  a digital  computer.  Simulation  signals  are  analyzed  by  the  practical 
techniques  presented  In  section  III  to  demonstrate  power-spectral  density  esti- 
mation, transfer  function  estimates,  coherence  and  confidence  calculations,  and 
frequency  averaging.  The  system  model  and  the  associated  numerical  equations  are 
developed  In  section  IV. 2.  Theoretical  analysis  results  are  derived  in  section 
IV. 3 and  the  simulation  and  analysis  results  are  presented  in  section  IV. 4. 

2 A MULTISENSOR  SYSTEM 

A suitable  numerical  example  system  must  be  sufficiently  generalized  to  ex- 
emplify the  significant  aspects  of  the  available  random  data  analysis  techniques. 
Simultaneously,  the  considered  system  must  be  simple  enough  to  be  amenable  to 
direct  theoretical  analysis  which  could  be  compared  to  the  simulated  results. 
Finally  we  desire  that  the  chosen  example  resemble  the  random  data  analysis  situ- 
ations which  we  encounter  In  practice. 

Most  generally  In  engineering  practice,  we  are  Interested  In  characterizing 
random  disturbances  and  the  responses  of  dynamic  systems  to  that  disturbance.  We 
are  constrained  to  practical  measurements  of  signals.  System  responses  and  dis- 
turbance phenomena  are  observed  with  Imperfect  sensors  which  contribute  additional 
random  components  (sensor  noise)  to  the  observed  or  measured  quantities.  For  ex- 
ample, the  measurements  and  analyses  of  random,  aerodynamlcally  Induced  pressure 
fluctuations  acting  upon  mechanical  structures  are  directly  affected  by  random 
signal  components  due  to  the  pressure  transducer,  the  Instrumentation  system,  and 
the  signal  processing,  as  well  as  the  desired  signal  directly  attributable  to 
actual  pressure  fluctuations. 


108 


AFWL-TR-76-193 


Two  distinct  cases  of  sensor  noise  data  corruption  occur.  In  one  case,  a 
sensor  Is  used  for  Instrumentation  purposes  only.  Sensor  noise  addltlvely  cor- 
rupts the  measured  quantity.  Another  case  of  Importance  occurs  with  sensor  Im- 
bedded as  a control  sensor  In  a dynamic  system.  System  dynamics  as  well  as 
sensor  noise  must  be  known  In  these  cases  to  properly  account  for  the  measured 
random  phenomena.  Examples  of  such  sensors  are  pressure  transducers  of  a hy- 
draulic actuator  or  the  tracking  sensor  of  a pointing  system. 

The  generalized  system  we  choose  to  Investigate  Includes  one  sensor  of  each 
type.  A simple  dynamic  system  described  by  a first-order  differential  equation 
Is  used.  Sensor  noise  Is  modeled  simply  as  white  noise  so  as  not  to  obscure  Im- 
portant aspects  of  random  data  analyses  with  a multitude  of  parameters.  Approxi- 
mation of  noise  spectra  by  white  noise  Is  valid  In  many  practical  cases  because 
the  noise  Is  "wlde-band"  in  comparison  with  system  dynamics.  Effects  of  "colored" 
sensor  noise  are  commented  upon  as  appropriate  In  the  theoretical  derivations. 

In  addition  to  control  loop  sensor  noise  the  dynamic  system  is  stochastically 
excited  by  an  exponentially  correlated  stochastic  disturbance  phenomena.  Such  a 
disturbance  model  Is  the  simplest  correlated  (colored)  disturbance  we  can  en- 
counter; yet  It  exemplifies  the  notions  of  correlatedness , bandwidth  and  non- 
constant spectra  which  we  typically  encounter,  though  often  in  a more  complicated 
manner.  In  practice. 

The  complete  system  model  used  for  our  theoretical  and  simulation  Investiga- 
tion Is  described  by  Its  transfer  function  block  diagram  given  in  figure  29.  The 
control  system  output  y Is  compared  with  the  desired  system  output  (reference  in- 
put) r by  the  control  system  sensor.  This  sensor  measures  the  control  error 
accurately  except  for  an  additive  sensor  noise  v.  The  noise  corrupted  sensor 
output  em  Is  the  only  system  error  signal  which  can  be  Instrumented  and  analyzed. 
The  error  signal  is  amplified  by  gain  K;  the  resulting  signal  drives  the  system 
plant  G.  The  system  plant  Is  also  driven  by  the  external  disturbance  d.  In  fact, 
the  purpose  of  the  control  system  Is  to  maintain  the  plant  output  y at  the  de- 
sired level  r (zero)  despite  perturbations  Induced  by  the  disturbance. 

The  Laplace  transform  complex  variable  s Is  used  In  the  representation  of  G 
and  the  stochastic  disturbance  model  (see  reference  6).  The  system  plant  Is  a 
simple  Integrator.  The  stochastic  disturbance  Is  modeled  as  the  output  of  a 


single-pole  low-pass  filter.  The  filter  has  bandwidth  A/2ir  (Hz)  and  unity  dc 
gain.  The  disturbance  model  excitation  and  the  sensor  noises  are  modeled  as 
Independent  normally  distributed  white  noise  stochastic  processes. 
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Figure  29.  Stochastic  Analysis  Example  Model. 

Model  Nomenclature 
K - control  system  loop  gain 
r - dynamic  system  reference  Input;  r(t)  = o 
y - control  system  output 
e - actual  system  error 

v - error  sensor  noise  (white  spectra,  gausslan  pdf) 
em  - measured  system  error 

d - exponentially  correlated  disturbance  phenomena 
w - disturbance  measurement  sensor  noise  (white  spectra,  gausslan  pdf) 
dm  - measured  system  disturbance 

n - disturbance  model  excitation  (white  spectra,  gausslan  pdf) 

X - disturbance  model  bandwidth  parameter  BW  (In  Hz)  = x/2tt 


In  absence  of  disturbance  phenomena,  the  control  system  of  figure  29  has 
zero  steady  state  error  and  output  responses.  Disturbances  excite  the  system  to 
give  non-zero  error  and  system  output.  RMS  characterizations  of  these  quantities 
are  a useful  figure  of  merit  Indicating  the  performance  of  the  system.  The 
smaller  the  RMS  of  these  quantities,  the  better  the  system  performance. 
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Equations  modeling  the  analysis  system  are  presented  In  the  remainder  of 
this  section.  In  section  IV. 3 theoretical  stochastic  relationships  among  the 
various  signals  are  developed.  Finally,  In  section  IV. 4 numerical  results  are 
presented.  Interpreted,  and  compared  with  the  theoretically  expected  results. 

The  disturbance  process  Is  modeled  by  a first  order  differential  equation. 

d = x(n-d)  (4-1) 

The  true  system  disturbance  Is  represented  by  d.  The  constant  A Is  a disturbance 
bandwidth  parameter  and  n Is  a white  noise  excitation  of  the  model.  Only  a noise 
corrupted  signal  Is  available  for  measurement.  The  noisy  disturbance  measurement 
Is  modeled  by  simply  adding  a white  noise  signal  component  w to  the  true  system 
disturbance.  Denoting  the  disturbance  measurement  by  dm  gives  the  following 
algebraic  relationship. 


The  control  system  Is  modeled  by  the  following  signal  relationships  evident  from 
figure  29 


y = Kem  + d = K(v-y)  + d 
and  the  algebraic  relationships 


(4-3) 


e = e + v 
m 


(4-4) 


and 


e = r-y  = -y  (r=o)  (4-5) 

Differential  equations  (4-1)  and  (4-3)  and  the  algebraic  equations  (4-2),  (4-4), 
and  (4-5)  specify  the  structure  of  the  system.  The  model  Is  completely  character- 
ized by  additionally  specifying  numerical  values  for  the  parameters  > and  K and 
by  specifying  the  statistics  of  the  white  noise  process  n,  w and  v.  These  white 
noise  processes  are  chosen  to  have  zero  means  and  spectral  densities  2,  4*^,  and 
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$vy,  respectively.  All  combinations  of  each  signal  having  either  a low  or  a high 
spectral  density  make  up  the  four  cases  considered.  Parameters  X and  K are  given 
values  2ir50  and  2irl0  respectively. 

The  equations  presented  here  are  not  suitable  for  numerical  simulation.  The 
differential  equations  must  be  approximated  by  finite-difference  equations. 

These  equations  are  obtained  by  approximating  derivatives  with  a finite  dlffer- 


x*^  Cx(t)  - x(t-At)] 


(4-6) 


Substituting  this  approximation  Into  differential  equations  (4-1)  and  (4-3)  leads 
to  the  finite  difference  equations 


dk  * dk-l  + AT  x(nk-l  " dk-l} 


(4-7) 


yk  - yk_!  + *T  K(uk_1  - yw)  ♦ AT  d^ 


(4-8) 


The  subscript  k denotes  the  successive  Increments  of  the  discrete  process  and  re- 
lates to  time  t = k aT  In  the  continuous  process. 


dk  ss  d(kAT) 


(4-9) 


The  discrete  data  algebraic  equations  are  simply: 


'V*k 


<em>k  ' ek  * uk 


(4-11) 


ek  ” -yk 
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Finally,  the  'white'  noise  sequences  nk>  vk,  and  wk  are  obtained  from  a random 
number  generator.  Samples  from  the  generator  have  variances  Qn/AT,QV/AT,  and 
C^/aT  respectively  and  are  Independent  of  each  other  and  from  sample  to  sample. 
Normalization  of  the  variances  by  aT  Is  required  for  an  accurate  approximation 
to  'continuous'  white  noise  which  has  Infinite  variance.  The  random  number  gen- 
erator described  In  reference  16  Is  used  In  the  simulation.  It  produces  Inde- 
pendent samples  having  normal  distributions.  Three  separate  starting  seeds  are 
used  to  produce  the  three  Independent  sequences. 

The  finite-difference  equations  used  to  produce  signal  sequences  for  stochas 
tic  signal  analysis  are  summarized  In  table  8.  These  equations  approximate  the 


Table  8. 

SUMMARY  EQUATIONS  FOR  THE  DISCRETE  TIME 
STOCHASTIC  SYSTEM  NUMERICAL  EXAMPLE 

DISTURBANCE  PROCESS  MODEL 

dk  • dk_,  ♦ »T  • x . (nk_,  - dk_,) 


CONTROL  SYSTEM  MODEL 

yk  »yk.,  YdT  • K(vk_,  -yk_,>  ♦ at  dk_, 


SYSTEM  ERROR  SENSOR  MODEL 

<en,>k  ' 9k  + vk  ’ - yk  * vk 


DISTURBANCE  SENSOR  MODEL 

ldA  ■ dk  + «k 


White  noise  sequences  nk>  vk>  and  wk  generated  by 
a gausslan  number  generator  to  have  variances  Qn/AT, 

Qv/aT,  and  Q^aT,  respectively. 

differential  equations  describing  a single-pole  control  system  excited  by  feed- 
back sensor  noise  and  a stochastic  disturbance.  The  approximation  Is  very  good 


in 


I 


1 
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at  low  frequencies;  It  suffers  from  aliasing  like  distortions  at  higher  frequen- 
cies. The  discrete  time  model  has  greater  phase  lag  and  less  magnitude  attenua- 
tion at  frequencies  approaching  one  half  the  discrete  equation  update  frequency. 
The  equations  given  In  Table  8 represent  a complete  mathematical  description  of 
the  process  under  Investigation.  In  section  IV. 3,  the  theoretical  analyses  of 
section  II  are  applied  to  determine  the  PSD's,  transfer  function,  and  coherence 
functions  we  expected  to  calculate  from  the  data.  These  theoretical  analyses  are 
possible  In  this  case  since  we  have  perfect  knowledge  of  the  system  and  the  sta- 
tistics of  the  component  stochastic  signals. 

3 THEORETICAL  COMPUTATIONS 

Given  the  mathematical  system  model  characterized  by  differential  equations 
(4-1)  and  (4-3)  along  with  algebraic  relationships  equations  (4-2),  (4-4),  and 
(4-5),  and  applying  to  linear  system  relationships  summarized  In  table  1 and 
equation  (2-8)  we  can  make  theoretical  calculations  of  signal  PSD  and  coherence 
functions  we  should  expect  to  calculate  by  random  data  analysis  procedures. 

These  theoretical  results  are  compared  with  the  signal  analysis  results  as  the 
latter  are  developed  in  section  IV. 4. 

Two  aspects  of  the  stochastic  analysis  model  example  are  Investigated. 

First,  the  disturbance  process  model  Is  studied;  PSD's  of  the  disturbance  sensor 
output  signal  dm  are  computed  for  several  signal -to-nolse  ratios  obtained  by  vary- 
ing the  amplitude  of  the  wideband  sensor  noise  signal  w.  Transfer  function  esti- 
mates for  the  disturbance  process  d/n  are  estimated  from  the  quantity  n (not  mea- 
surable In  practice)  and  the  measured  disturbance  dm>  Transfer  function  methods 
for  several  RMS  levels  of  w are  compared.  Secondly,  transfer  function  estimates 
as  well  as  PSD  and  coherence  calculations  are  made  from  the  measurement  of  system 
error  em  and  the  measured  system  disturbance  dm.  These  quantities  are  also  com- 
puted for  a variety  of  cases  In  which  the  signal  to  noise  ratio  of  each  sensor  Is 
varied  over  a lOOdB  range. 

The  basic  analysis  tool  Is  the  transfer  function  H(2irf).  Two  sybsystems  of 

the  numeric  example  are  analyzed.  We  first  consider  the  simple  process  relating 

measured  disturbance  dm  to  the  white  noise  excitation  n.  We  realize  that  a slm- 

m 

liar  analysis  Is  not  possible  with  a practical  system  since  no  signal  correspond- 
ing to  n Is  physically  available.  However,  we  proceed  with  analysis  of  this 
simple  process  as  It  shall  lend  Insight  to  the  use  and  Interpretation  of  our  pro- 
cedures. Next  we  shall  consider  the  process  relating  system  measure  error  em  with 
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the  measured  disturbance  d . These  functional  relationships  are  depicted  ex- 

m 


pllcltly  In  block  diagrams  of  figure  30. 


X*  2t50 


Signals  Available 
for  Analysis 
n 
d 

dm 


a.  Disturbance  Process  Model 


Signals  Available 
for  Analysis 

dm 

em 


K = 2t1  0 


Figure  30.  Theoretical  Process  Transfer  Functions. 
The  transfer  function  from  Input  n to  output  d Is  seen  to  be 

«»>  ■ rW 

It  follows  from  equation  (2-80)  that 

*dd(f)  - |H(2Wf)l2*nn(f) 

•[l  +(¥)2l"1*nn(f) 


(4-13) 


(4-14) 
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The  additive  noise  w adds  spectral  component  « to  the  measured  disturbance 

WW 

signal.  Hence,  the  relationship 


d = d + w 
m 


(4-15) 


Implies 


d d 
m m 


*dd  + 


ww 


(4-16) 


because  the  processes  n and  w have  been  assumed  (and  constructed)  to  be  indepen- 
dent. A useful  representation  of  4dmdm  equivalent  to  that  of  equation  (4-16)  is 
obtained  by  factoring 


i 


3 *dd 
m m 


(4-17) 


The  second  term  within  the  brackets  Is  a noise-to-slgnal  ratio,  which  we  denote 
notatlonally  as  N/S.  Observe  that  N/S  is  a function  of  frequency  since,  In  gen- 
eral, both  * and  «dd  are  functions  of  f.  Over  those  frequency  bands  In  which 
N/S  « 1,  the  measurement  spectrum  closely  approximates  the  spectrum  of  the  de- 
sired signal  d.  The  noise  Is  Inconsequential  In  these  cases.  Over  the  remaining 
frequency  bands,  the  noise  term  contributes  significantly  to  the  measured  spectrum. 
Substituting  equation  (4-14)  Into  equation  (4-17)  gives  the  result 


\Vf)  ■ «J”<^>I2 


i ♦ 


ww 


|H(2-rrf ) |2  $ 


nn  J 


(4-18) 


This  PSD  Is  plotted  In  figure  31  for  several  ratios  of  ®ww/®nn*  In  the  examples 
considered,  *nn  and  are  constants  because  they  represent  spectral  densities 
of  'white'  noise  processes.  Thus,  a low  slgnal-to-noise  condition,  as  previously 
defined,  exists  whenever 


|H(2irf)|2 


ww 


"nn 


(4-19) 
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Figure  31.  Disturbance  Measurement  Power  Spectral  Densities. 

In  addition  to  PSD's,  we  are  also  Interested  In  transfer  function  estimates  and 
the  coherence  function.  These  functions  are  computed  by  use  of  the  equations  de- 
veloped In  section  1 1. 9.  The  measured  disturbance  signal  Is  taken  as  the  trans- 
fer function  output.  The  Input  Is  the  white  noise  excitation  n (see  figure  30a). 

TRANSFER  FUNCTIONS 

METHOD  1 

H(w)  - *d  n (w/2ir)/*nn 
m 

= X+3w  (4-20) 

Method  1 yields  an  Ideal  estimate.  The  output  noise,  as  long  as  it  Is  Incoher- 
ent with  the  system  excitation,  has  no  effect  on  the  transfer  function  estimate. 
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METHOD  2 


I H (w ) | 


r*d  d (w/2ir)/$  ] i/2 

*■  mm  J 

1/2  i1  * f 1 + 


(r)1 


* jl/2 

ww 

$ 

vnn 


(4-21) 


Thus  Method  2 estimates  the  correct  transfer  function  in  this  case  only  when  the 
output  signal  measurement  as  well  as  the  Input  measurement  is  noise  free.  The 
Method  2 estimate  is  plotted  in  figure  32. 


Figure  32.  Disturbance  Process,  Method  2 
Transfer  Function  Estimates. 


COHERENCE  FUNCTION 

The  coherence  function  y2(f)  Is  calculated  by  using  equation  (2-109) 
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YVf) 

m 


•dn(f)  2 
m 

\d(f)  *nn 
mm 


$ 

ww 


l-i 


\ - 2ir50  (4-22) 

We  see  that  If  $ « $ that  y2  Is  approximately  unity  at  frequencies  less 

ww  nn 

than  50  Hz.  At  higher  frequencies,  y2  becomes  smaller.  The  theoretical  coher- 
ence function  Is  plotted  in  figures  39  and  40  of  section  IV. 4 in  comparison  with 
a coherence  function  calculated  from  the  actual  data.  Observe  that  the  coher- 
ence function  is  near  unity  over  those  frequency  bands  in  which  both  transfer 
function  estimates  are  accurate.  Also,  in  the  same  frequency  bands  the  coherent 
output  power  dominates  the  sensor  noise. 


Identical  procedures  are  employed  to  calculate  the  power  spectral  density 
<J»em  em-  In  this  linear  system,  the  error  em  is  made  up  of  two  Independent  com- 
ponents - one  due  to  the  disturbance  d and  the  other  due  to  the  sensor  noise  v. 
The  appropriate  transfer  functions  are 


em  a 
d“  = 


Mw) 


1/K 

1 + J'w/K 


(4-23) 


and 


v31  - Hi(w)  = r^jw/R'  = jwHl  (w) 


(4-24) 


By  superposition  and  from  the  linear  system  theorem  we  have 


Vm  = IM2»OI2  «dd+|H2(2.f)l2  ♦ 


vv 


(4-25) 


Notice  that  In  the  latter  case,  the  precise  manner  by  which  the  error  sensor 
noise  corrupts  the  error  measurement  is  Influenced  by  the  control  system  dynamics. 
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This  is  a general  result  true  for  any  sensor  embedded  as  a feedback  sensor  in  a 

control  system.  Plots  of  the  theoretical  e are  presented  in  figure  33. 

m m 


Figure  33.  Control  Error  Measurement  PSD's. 


The  measured  signals  dm  and  em  may  be  considered  as  the  Input  and  output, 
respectively,  of  a linear  system.  These  signals,  each  corrupted  by  additive 
noise  components,  represent  the  general  signal  analysis  transfer  function  estima- 
tion problem  encountered  in  practice.  Theoretical  transfer  function  estimates 
based  on  the  model  transfer  function  and  the  noise  characterizations  are 

IDEAL 


= K 1 +1  jw/K 


(4-26) 
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METHOD  1 


H(w)  = 


*e  d 
m m 

fdd 
m m 


H(wJ* 


dd 


'<L<L 
m m 


(4-27) 


= H(w) 


(4-28) 


The  last  step  follows  from  equation  (4-17).  We  again  note  the  Importance  of  a 
noise- to-slgnal  ratio,  N/S.  In  this  case  the  noise  Is  characterized  by  * 

WW 

while  $dd  characterizes  the  signal.  Over  low  N/S  frequency  bands  the  transfer 
function  estimate  Is  accurate.  Conversely,  over  high  N/S  frequency  bands  the 
estimate  Is  poor.  Observe  that  the  estimated  transfer  function  magnitude  is 
less  than  the  true  magnitude.  Also  note  that  the  estimated  transfer  function 
Is  a real  scalar  multiple  of  the  actual  transfer  function.  Consequently,  the 
phase  component  of  the  estimate  remains  accurate  despite  the  noise-to-slgnal 
level . 

METHOD  2 


|H(w)| 


emem 

*d  d 
m m 


1 1/2 


14W*  ^ 

*dd 


-I/2 


H(w) 


1 + 


WW 

*dd 


(4-29) 


Only  the  transfer  function  magnitude  can  be  estimated  by  Method  2.  As  for  the 
Method  1 estimate,  this  estimate  magnitude  differs  from  the  true  transfer  func- 
tion magnitude  by  a scalar  multiple.  In  this  case,  the  multiple  may  be  either 
less  or  greater  than  unity.  The  numerator  terms  represent  an  output  signal 
noise-to-slgnal  ratio.  The  w2  factor  appears  In  this  term  to  compensate  for  the 
fact  that  the  signals  characterized  by  4>vy  and  $dd  occur  In  different  locations 
of  the  control  system;  to  properly  compare  them,  they  must  be  transferred  to  a 
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common  comparison  point.  The  output  (see  figure  30a)  Is  chosen  as  the  common 
comparison  poln?.  Observe  that  the  relative  transfer  function  from  d or  v to 
the  output  em  Is  simply  1/s  corresponding  to  H2  ■ 1/w2.  This  factor  Is  the 
appropriate  multiplier  of  0^  so  that  It  can  be  compared,  on  a nolse-to-slgnal 
basis,  with  *w.  Similarly,  the  denominator  factor  also  has  a N/S  term.  An 
accurate  transfer  function  estimate  Is  obtained  whenever  the  N/S  ratios  are 
much  less  than  unity.  For  larger  N/S  ratios,  the  estimated  transfer  function 
deviates  from  the  actual  transfer  function  (except  for  the  pathological  case 
wherein  the  numerator  and  denominator  N/S  ratios  vary  so  as  to  maintain  a unity 
multiplying  error  factor). 

COHERENCE  FUNCTION 


Y*d  e 
m m 


Vm  Vn, 


-I 


(4-30) 


I 


■ 

i 


Observe  that  the  coherence  function  Is  made  up  of  N/S  factor  terms.  As  either 
transfer  function  estimate  (Method  i or  2)  deviates  from  the  actual  transfer  func- 
tion, y2  reduces  In  value  from  unity.  Also  note  that  y is  simply  the  ratio  of 
the  Method  1 to  the  Method  2 transfer  function  amplitudes. 

The  transfer  function  estimates  and  the  coherence  function  theoretical 
values  are  plotted  In  comparison  with  transfer  functions  and  coherence  functions 
calculated  from  the  simulated  model  data  In  section  IV. 4.  We  note  the  general  re- 
sult that  typically,  as  the  system  noise  Increases,  these  functions  diverge  from 
their  desired  values. 

4.  SIMULATION  AND  ANALYSIS  RESULTS 

Results  of  applying  the  numerical  analysis  procedures  developed  In  sections 
II  and  III  to  the  simulated  control  system  using  the  SIGANAL  software  are  pre- 
sented. These  results  are  compared  with  theoretical  values  to  demonstrate  repre- 
sentative accuracies  and  Interpretative  procedures  useful  In  analyzing  PSD  and 
transfer  function  estimate  data. 
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For  the  hypotehtlcal  system  used  in  this  simulation,  two  comparisons  of  the 
results  are  possible  and  equally  instructive.  The  analysis  results  are  compared 
to  Ideal  values.  Transfer  function  estimates,  for  example,  are  compared  to  the 
actual  system  transfer  functions.  Also,  analysis  results  are  compared  to  those 
theoretically  expected  results  as  derived  In  section  IV. 3.  This  later  compari- 
son Is  only  possible  In  this  simulated  example,  since  the  theoretical  calcula- 
tions require  precise  knowledge  of  system  transfer  functions  and  measurement 
noise  properties  which  are  not  known  In  practice. 

The  system  was  simulated  by  Iterating  the  modeling  equations  In  table  8. 

The  numeric  values  obtained  for  the  disturbance  process  nk,  d^,  (dm)k.  wk>  and 
the  control  system  v^,  yk  and  (em)k  were  written  In  blocks  on  an  engineering 
unit  (EU)  format  data  tape  for  subsequent  analysis  by  the  SIGANAL  software  (ref- 
erence 1). 

The  simulated  signals  were  obtained  with  a finite  difference  update  rate  of 
1 k Hz  corresponding  to  a simulation  of  time  Increment  aT  = .001  second.  The 
equations  were  Iterated  to  give  8192  signal  values  simulating  T = 8.2  seconds  of 
system  response. 

The  simulated  signals  stored  on  EU  tape  are  analyzed  by  the  SIGANAL  software. 
The  SIGANAL  summary  list  and  NAMELIST  Input  cards  used  for  the  analysis  are 
listed  In  table  9.  Excerpts  of  these  analysis  results  are  presented  and  dis- 
cussed In  the  remainder  of  this  section. 

a.  Disturbance  Process  PSD's 

The  results  of  the  disturbance  process  analysis  presented  in  figures  34-40 

are  first  reviewed.  The  disturbance  process  excitation  signal  measured  PSD  Is 
plotted  In  figures  34a  and  b.  These  measured  values  differ  somewhat  from  the 
actual  value  *nn  * 2.0  also  plotted  because  of  the  statistical  uncertainty  of 
PSD  estimates  made  from  a finite  duration  observation  of  the  signal  as  discussed 
In  section  III. 4.  Confidence  bounds  for  95%  confidence  are  also  plotted  and  were 
derived  by  use  of  the  methods  Illustrated  in  Example  3-2.  The  PSD  estimate  In 
figure  34a  Is  obtained  from  a minimum  of  eight  frequency  averages  while  the  PSD 
In  figure  34b,  which  varies  more  from  the  theoretical  value,  is  obtained  with  a 
minimum  of  only  four  averages.  The  95%  confidence  bound  calculation  results  are 
summarized  In  table  10  for  each  of  these  cases.  Case  A employs  at  least  four 
times  the  number  of  raw  PSD  estimates  per  average  than  Case  B.  The  spread  of  the 
SIGANAL  PSD  estimates  scale  with  the  confidence  bound  calculations.  Comparisons 
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Figure  35.  True  Disturbance  PSD. 


Disturbance  Process  Cross-Spectral  Densities. 


Disturbance  Process  CSD  Phase 


39.  Disturbance  Process  Transfer  Function  Estimates,  4>  = .02. 


39  (Continued).  Disturbance  Process  Transfer  Function  Estimates, 


Figure  40.  Disturbance  Process  Transfer  Function  Estimates, 
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Table  10 


STOCHASTIC  SYSTEM  NUMERICAL  EXAMPLE 
95%  CONFIDENCE  BOUNDS  FOR  PSD  ESTIMATES 


Frequency 

Band 

(Hz) 

Averaged 

Resolution 

Bandwidth 

(Hz) 

Degrees  of 
Freedom 
N = 2(BW/aF) 

Confidence 

Multipliers 

M M 

i V 

Spread 

(magnitude) 

Case 

0-20 

1 

16 

.4 

1.8 

4.1 

Case 

A 

20-100 

4 

65 

.7 

1.4 

2.0 

A 

100-250 

10 

164 

.8 

1.2 

1.6 

250-500 

25 

410 

.9 

1.1 

1.3 

Case 

0-20 

.25 

4 

.1 

2.8 

23 

Case 

B 

20-100 

.50 

8 

.3 

2.2 

8.1 

B 

100-250 

2.5 

41 

.6 

1.6 

2.4 

of  figures  34a  and  34b  show  dramatically  the  advantage  of  adequate  frequency 
averaging  to  obtain  smooth  PSD  estimates. 

Similarly,  the  true  disturbance  process,  d,  output  PSD  estiamte  4^  is 
plotted  In  figures  35a  and  b along  with  the  theoretical  value  and  the  95%  con- 
fidence bounds.  Observe  that  the  calculated  values  closely  follow  the  trend  of 
the  theoretical  value  and  the  95%  confidence  bounds.  Further  note  that  the  cal- 
culated PSD  varies  about  the  theoretical  value  with  amplitudes  that  scale  with 
the  confidence  bounds.  This  result  again  demonstrates  the  advantage  of  adequate 
frequency  averaging;  the  actual  PSD  is  readily  more  apparent  from  the  analysis 
calculations  In  figure  35a  than  the  similar  calculation  with  less  averaging 
shown  in  figure  35b. 

A second  theoretical  PSD  curve  has  been  plotted  in  figure  35a  to  illustrate 
the  distortions  from  the  ideal  introduced  by  the  finite-difference  equations  used 
to  model  the  process.  The  distortion  causes  less  amplitude  attenuation  than  the 
Ideal  filter  at  frequencies  approaching  1/2aT  as  explained  in  section  IV. 2.  The 
correspondence  of  the  SIGANAL  PSD  estimate  with  the  finite-difference  theoretical 
PSD  shows  that  the  SIGANAL  calculations  estimate  accurately  the  actual  signal  PSD. 
The  finite-difference  PSD  Is  expressed  by  the  equation 

* 3 2(UT)2  | [cos  ( 2irf At ) - (1  - XAT)]2  + [si n ( 2irf aT ) ]2 (4-31) 
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which  Is  obtained  from  the  finite  difference  equation,  equation  (4-7),  by  methods 
described  in  references  6,  10,  11,  17,  and  18. 

The  disturbance  process  Is  a low-pass  process.  At  frequencies  below  10  Hz 
the  white  noise  excitation  Is  passed  by  the  filter  without  attenuation  (compare 
figures  34a  and  34b).  At  frequencies  above  10  Hz,  the  process  attenuates  the  In- 
put by  a factor  that  Increases  with  frequency  as  described  by  equation  (4-14). 

This  example  demonstrates  the  filtering  action  of  linear  systems. 

b.  Disturbance  Process  CSD's 

Given  the  disturbance  process  excitation  n and  the  output  y,  cross-spectral 
densities  are  computed.  The  CSD  amplitude  is  simply  the  amplitude  of  the  Input 
spectra,  in  this  case  a constant  4>nn  = 2.0,  multiplied  by  the  amplitude  of  the 
disturbance  process  transfer  function.  The  SIGANAL  CSD's  are  presented  In  figures 
36a  and  b where  different  frequency  averaging  schedules  have  been  used  as  In  the 
previous  figures  and  as  listed  In  table  10.  The  theoretically  expected  CSD  ampli- 
tude Is  followed  by  the  SIGANAL  calculations  with  only  statistical  variations  due 
to  the  finite  data  observation  Interval.  Like  the  PSD,  the  CSD  estimate  improves 
with  more  averaging.  Exact  confidence  bounds  for  CSD's  have  not  been  calculated, 
but  figure  36  suggests  that  they  are  approximately  the  same  as  PSD  bounds. 

The  disturbance  process  Input  output  cross-spectral  density  phase  is 
plotted  In  figure  37.  The  same  phase  curve  Is  obtained  independently  of  the  fre- 
quency averaging  used.  The  SIGANAL  CSD  phase  calculation  deviates  somewhat  from 
the  Ideal  phase  curve  but  follows  exactly  the  phase  of  the  finite-difference  nu- 
meric model  As  in  the  rase  of  thp  disturbance  precede  output  PSD,  differences 
between  the  Ideal  and  the  calculated  phase  values  stem  from  the  distortions  of  the 
finite-difference  simulation  equations  rather  than  from  Inaccuracies  of  the  CSD 
phase  calculation.  For  reference,  the  Ideal  phase  + Is  obtained  from  the  process 
transfer  function  equation  (4-13)  as 


♦ideal  = 'tan'1  (¥■) 


(4-32) 


while  the  phase  associated  with  the  finite-difference  model  Is 


*F-D  =“tan_1 


s 1 n ( 2-rf  aT  ) 

cos(2rfAT)-  (V-  XaT) 


(4-33) 
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The  CSD  phase  calculation  does  not  exhibit  the  statistical  variation  we 
saw  In  the  PSD  and  CSD  magnitude  estimates.  This  result  Is  typically  valid  when- 
ever the  signals  used  to  generate  the  CSD  estimate  represent  the  actual  Input  and 
outputs  of  a linear  system  and  that  further,  neither  measurement  Is  corrupted  by 
noise.  Examples  subsequently  presented  (e.g.,  figure  40d)  show  that  the  CSD 
phase  signal  does  exhibit  variations  from  the  true  linear  system  phase  character- 
istic when  one,  or  both,  of  the  signals  become  corrupted  by  an  additive  noise 
signal . 

c.  Disturbance  Process  Transfer  Function  Estimates 

Transfer  function  estimates  of  the  disturbance  process  are  obtained  from 
the  Input  and  output  PSD*s  and  their  CSD.  The  Method  1 transfer  function  ampli- 
tude Is  the  CSD  magnitude  normalized  by  the  Input  excitation  PSD,  $nn-  The  re- 
sulting Method  1 transfer  function  estimate  is  plotted  in  figure  38a.  The  Method 
2 estimate  obtained  by  taking  the  square  root  of  the  ratio  of  the  output  PSD  to 
the  Input  PSD  as  In  equation  (4-21)  matches  exactly  the  Method  1 estimate  since 
there  are  no  extraneous  noise  sources  present.  Both  estimates  deviate  from  the 
Ideal  transfer  function,  also  plotted,  but  match  exactly  the  actually  simulated 
transfer  function.  Figure  38a  shows  that  precise  transfer  function  estimates  can 
be  calculated  by  either  Method  1 or  Method  2 techniques. 

The  transfer  function  phase  is  simply  the  CSD  phase  in  Method  1.  No  tech- 
nique for  calculating  the  transfer  function  phase  is  presented  in  the  Method  2 
analysis.  Those  familiar  with  linear  systems  theory  might  suggest  that  a phase 
estimate  could  be  obtained  by  differencing  the  Fourier  transform  phases  of  the 
system  output  and  Input.  The  phase  estimate  obtained  in  this  manner  is  plotted 
In  figure  38b  along  with  the  ideal  phase.  The  estimated  phase  Is  extremely  noisy, 
especially  at  higher  frequencies,  and  Is  attributable  to  the  phase  uncertainty  of 
Fourier  transforms  of  stochastic  processes,  even  when  there  are  no  additive  noise 
signals.  Because  of  their  high  uncertainty,  phase  estimates  obtained  In  this 
manner  are  seldom  used  in  practice.  We  shall  restrict  our  attention  henceforth  to 
CSD  phase  estimates. 

The  coherence  function  Is  not  plotted  for  this  no  noise  case.  With  no 

noise,  the  nolse-to-slgnal  ratio  terms  vanish  [e.g.,  * - 0 In  equation  (4-22)] 

ww 

and  the  Ideal  coherence  function  Is  unity  over  the  entire  frequency  band.  The  co- 
herence calculated  by  the  SIGANAL  routines  Is  not  precisely  unity  because  of  nu- 
merical roundoff.  However,  the  calculated  coherence  values  are  In  the  range 
(.998,  1]  which  Indicates  extremely  high  coherence  between  the  process  Input  and 
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output  signals.  We  expect  high  coherence  because  the  output  signal  Is  totally 
generated  by  application  of  the  Input  signal  to  the  disturbance  process. 

d.  Noise  Corrupted  Disturbance  Process 

The  previous  example  dealt  with  a noise  free  measurement  of  the  process 
Input  and  output  signals.  In  the  following,  only  the  noise  corrupted  process  out- 
put signal  dm  is  available  for  analysis.  The  additive  noise  signal  w adds  biases 
and  variations  to  signal  PSD  and  transfer  function  estimates.  Transfer  function 
estimates  by  Methods  1 and  2,  coherence  calculations  and  the  CSD  phase  are  plotted 
In  figures  39a  - d.  In  this  case,  * equals  .02  corresponding  to  a noise  to  sig- 

WW 

nal  (N/S)  ratio  of  1/100  in  the  low  frequency  band  (<  50  Hz).  Note  that  the  trans- 
fer function  magnitude  estimates  display  a more  erratic  character  than  the  noise 
free  estimates  In  figure  38.  The  presence  of  corruptive  noise  Is  Indicated  by  the 
coherence  function  which  Is  biased  from  unity  and  which  diverges  further  from 
unity  at  frequencies  beyond  50  Hz.  This  behavior  Is  predicted  by  a N/S  analysis. 

The  constant  amplitude  noise  ® * .02  becomes  a larger  fraction  of  the  signal 

PSD  which  Is  increasingly  attenuated  at  higher  frequencies. 

The  degrading  effect  of  higher  N/S  ratios  Is  demonstrated  In  figure  40.  A 
low  frequency  band  N/S  ratio  of  1 is  simulated  by  Increasing  ® to  2.0.  Both 
Method  1 and  Method  2 transfer  function  amplitude  estimates  and  the  CSD  phase  be- 
come appreciably  erratic.  The  Method  2 transfer  function  estimate  has  also  be- 
come biased  from  the  true  value.  This  bias  arises  from  the  output  measurement 
noise  which  adds  a bias  to  the  true  process  output  PSD  as  shown  In  equation  (4-16). 
The  biased  output  PSD  biases  the  transfer  function  estimate.  The  presence  of 
additive  noise  components  corresponding  to  large  N/S  ratios  is  indicated  Dy  the 
coherence  function,  figure  40c  which  has  dropped  to  at  most  .5.  As  shown  In  sec- 
tion III,  a coherence  of  .5  shows  that  only  half  of  the  output  signal  power  is 
attributable  to  a linear  system  excited  by'the  input  signal.  Thus,  as  Inspection 
of  figure  30a  shows,  $ contributes  to  the  measured  output  power;  ® contributes 
one-half  the  total  output  power  In  this  case. 

e.  Control  System  PSD's 

The  control  system  analysis  results  are  presented  in  figures  41-44.  A sum- 
mary of  the  theoretical  transfer  function  estimates,  which  are  well -matched  by  the 
SIGANAL  results,  are  plotted  In  figures  45-47. 
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Figure  41.  Control  System  Error  PSD's. 


Figure  42.  Control  System  Transfer  Function  Estimates  $ = 2E-5  * = *02. 


gure  43.  Control  System  Transfer  Function  Estimates  * = 2E-5 


(Continued).  Control  System  Transfer  Function  Estimates 
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Figure  47.  Coherence  Function  Summary. 


The  power  spectrum  of  the  control  system  error  depends  upon  the  error  per- 
turbation of  the  system  disturbance  input  d and  the  component  due  to  the  sensor 
error.  The  contribution  to  the  error  signal  PSD  of  each  of  these  sources  Is  char- 
acterized by  equation  (4-25)  developed  in  section  IV. 3.  These  theoretical  re- 
sults are  plotted  In  figure  33.  The  PSD  estimate  obtained  by  analyzing  the  digi- 
tally simulated  error  signal  with  the  SIGANAL  routines  Is  plotted  in  figure  41. 

The  theoretical  PSD  estimates  of  figure  33  are  also  plotted  to  show  the  excellent 
agreement  of  the  calculated  and  theoretical  results.  The  error  PSD  $emem  f°r  l°w 
sensor  noise  (*yy  = .00002)  Is  plotted  in  figure  41a  and  for  high  sensor  noise 
(#yy  a • 002)  in  figure  41b.  The  low  frequency  behavior  of  $emem  unalterated  by 
these  variations  In  4yy  while  the  high  frequency  PSD  amplitude  varies  directly 
with  the  changes  In  «yy,  the  error  sensor  noise.  This  observation  suggests  the 
hypothesis  that  the  low  frequency  error  contribution  Is  dominated  by  the  dis- 
turbance process  and  that  the  high  frequency  behavior  is  dominated  by  the  error 
sensor  noise.  The  hypothesis  Is  correct  for  this  simulated  control  system  and 
the  ranges  of  «yy  considered.  Moreover,  in  the  $yy  = .00002  case,  the  error  PSD 
shows  that  the  system  error  energy  Is  concentrated  in  the  low  frequencies  (<  20  Hz 


I 


146 


AFWL-TR-76-193 


or  so).  When  $vv  Is  Increased  to  .002,  the  control  system  error  becomes  concen- 
trated in  the  high  frequency  band  (>  20  Hz).  Comparisons  of  figures  41a  and  41b 
show  that  the  RMS  of  the  error  Is  much  larger  in  the  latter,  since,  by  equation 
(2-111),  the  RMS  Is  given  by  the  square  root  of  the  integral  of  4>emem  over  the 
entire  frequency  band,  and  since  the  PSD  in  figure  41b  is  greater  in  magnitude 
than  the  PSD  of  figure  41a. 

f.  Control  System  Transfer  Function  Estimates 

Estimates  of  the  closed  loop  transfer  function  from  the  disturbance  exci- 
tation (Input)  to  the  system  error  (output)  are  obtained  by  appropriately  pro- 
cessing measurements  of  the  system  disturbance  and  the  system  error.  These  esti- 
mates are  obtained  by  using  Methods  1 and  2 transfer  function  formulae  as  had  been 
used  to  obtain  disturbance  process  transfer  function  estimates.  The  control  pro- 
cess transfer  function  estimate  problem  adds  several  additional  features  not  ex- 
hibited In  the  disturbance  process  example.  They  are: 

• The  system  input  signal  measurement,  as  well  as  the  output,  is  cor- 
rupted by  additive  measurement  noise. 

• The  system  Input  signal  PSD  varies  with  frequency  (the  disturbance  is 
a "colored"  or  correlated  process). 

• The  output  measurement  noise  Is  Influenced  by  the  control  system  In 
which  the  noisy  error  (output)  sensor  is  embedded. 

The  transfer  function  estimates  are  presented  In  figures  42-44  for  three 
cases  of  noise  level.  In  figure  42  the  case  of  low  measurement  noise  is  examined. 
The  high  input  measurement  noise  case  (with  low  output  noise)  is  examined  In  fig- 
ure 43.  High  output  noise  but  low  input  noise  is  considered  in  figure  44.  Method 
1 transfer  function  SIGANAL  results  are  plotted  In  part  a of  each  figure,  Method  2 
results  are  plotted  in  part  b and  the  coherence  function  results  are  plotted  In 
part  c.  Also  plotted  in  each  figure  are  the  theoretical  results  that  we  would  ex- 
pect to  obtain  given  perfect  knowledge  of  the  system  transfer  functions,  the  true 
disturbance  PSD  and  the  measurement  noise  PSDs.  These  theoretical  results  are 
plots  of  equations  (4-28),  (4-29),  and  (4-30)  derived  in  section  IV. 3.  These 
theoretical  results  are  summarized  in  figures  45-47.  Also  plotted  with  the  trans- 
fer function  estimates  Is  the  Ideal  (or  actual)  transfer  function,  which  typically 
differs  from  the  theoretical  transfer  function.  Differences  between  the  ideal  and 
the  theoretical  transfer  functions  are  distortions  Introduced  into  the  theoretical 
transfer  function  estimate  by  the  noise  signals  Inherent  In  the  signal  measure- 
ments. These  distortions  are  systematic  errors  and  cannot  be  removed  from  the 
transfer  function  estimation  process. 
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Except  for  the  case  of  the  high  error  sensor  noise  Method  1 transfer  func- 
tion estimate  in  figure  44a,  the  SIGANAL  transfer  function  results,  obtained  by 
processing  the  simulated  signals,  agree  very  closely  with  the  theoretically  ex- 
pected results.  Thus,  the  SIGANAL  processing  is  shown  to  provide,  in  a practical 
data  analysis  situation,  accurate  results.  The  erratic  variations  of  the  computed 
estimates  about  the  theoretical  curves  are  variations  within  the  statistical  con- 
fidence as  limited  by  the  finite  duration  data  interval  processed.  These  varia- 
tions would  Increase  if  less  frequency  averaging  were  employed  to  smooth  the  raw 
estimates.  Conversely,  the  variations  could  be  reduced  by  frequency  averaging 
more  adjacent  estimates  to  Increase  the  number  of  degrees  of  freedom.  The  esti- 
mates could  also  be  improved  by  processing  a longer  data  interval.  The  Case  A 
frequency  averaging  schedule  listed  in  table  10  is  employed  to  obtain  the  esti- 
mates plotted  in  figures  41-44. 

The  anomalous  variation  of  SIGANAL  and  theoretical  results  in  figure  44a 
stems  from  a low  statistical  confidence  of  the  SIGANAL  result  in  the  frequency  band 
of  30  Hz  and  above.  Deviation  of  statistical  confidence  bounds  for  transfer  func- 
tion estimates  is  extremely  complex.  Is  data  (transfer  function)  dependent,  and 
often  cannot  be  practically  and  accurately  calculated  in  a given  problem.  Nonethe- 
less, we  argue  that  the  variation  of  results  exists  because  of  low  statistical  con- 
fidence caused  by  the  relatively  high  noise  level  of  the  error  signal  in  the  high 
frequency  band.  Reference  figure  41b.  The  cross-spectral  density  between  the 
Input  and  output  (error)  signals  In  this  high  frequency  band,  while  theoretically 
having  near  zero  magnitude,  computationally  has  components  large  compared  with  with 
the  true  CSD,  unless  enough  averaging  can  be  performed  to  improve  the  computational 
average.  Arguments  as  to  the  source  of  this  variation  are  academic  in  the  practi- 
cal data  analysis  case  where  we  must  ignore  or  disregard  transfer  function  esti- 
mates over  frequency  bands  in  which  the  coherence  is  low.  In  this  case,  the  co- 
herence is  low  above  5 Hz  and  we  must  ignore  the  transfer  function  estimates. 

The  coherence  function  Indicates  the  relative  level  of  noise  present  in 
either  the  Input  or  output  measurement.  The  coherence  function  decreases  as  the 
noise  level  increases.  Since  Increased  noise  levels  bias  and  corrupt  our  transfer 
function  estimates,  we  must  discount  these  estimates  in  high  noise  situations, 
which  we  can  determine  from  the  coherence  function.  As  a rule  of  thumb,  we  should 
not  trust  estimates  made  in  frequency  bands  for  which  the  coherence  is  0.5  or  less. 
Theoretically , output  signal  noise  does  not  effect  the  Method  1 transfer  function 
estimate.  Unfortunately,  there  is  no  technique  by  which  we  can  take  advantage  of 
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this  fact  since  we  cannot  determine  whether  the  input  or  the  output  signal  noise 
is  causing  low  coherence. 

Figures  45  and  46  show  that  the  Method  1 transfer  function  estimate  is  the 
better  estimate.  This  result  Is  generally  true  since  Method  1 estimates  are  sen- 
sitive to  only  input  signal  measurement  noise  while  Method  2 estimates  are  sensi- 
tive to  both  input  and  output  measurement  noises.  Examination  of  the  coherence 
functions  plotted  in  figure  47  shows  that  the  best  transfer  function  estimates  are 
obtained  at  low  frequencies.  At  higher  frequencies,  the  N/S  increase  thus  lower- 
ing the  coherence  and  corrupting  the  transfer  function  estimates. 

g.  Miscellaneous  Topics 

The  simulated  signal  analysis  example  provides  a mechanism  for  additional 
insight  and  analysis  that  might  be  applied.  The  character  of  stochastic  process 
signals  Is  illustrated  in  figures  48a  and  48b  in  which  the  time  histories  of  two 
simulated  signals  are  plotted.  The  control  system  output  signal  y Is  plotted  in 
figure  48a.  This  signal  Is  limited  in  frequency  content  to  10  Hz  by  the  band- 
width of  the  control  system.  Notice  the  slow  temporal  variations  of  the  signal. 
The  true  disturbance  signal  plotted  in  figure  48b  has  a noticeably  higher  fre- 
quency content.  This  signal  Is  the  output  of  a 50  Hz  process  which  does  contain 
higher  frequency  components.  A review  of  the  signal  time  histories,  as  these 
plots  provide,  should  be  performed  In  any  analysis  to  determine  whether  the  signal 
behavior  Is  random,  as  the  figure  48  plots  appear,  or  whether  there  are  determin- 
istic components  characterized  by  trends,  discontinuities,  or  low  frequency  varia- 
tions . 

Another  useful  analysis  is  a probability  density  function  calculation  of 
the  measured  signal  values.  The  probability  density  function  Indicates  the  rela- 
tive frequency  with  which  signal  values  occur.  The  pdf  also  characterizes  signal 
mean  values  and  the  dispersion  or  variance  of  the  signal  values  about  the  mean. 
Finally,  the  pdf  allows  the  analyst  to  determine  If  the  density  function  Is  gaus- 
slan,  symmetric,  or  skewed  and  whether  or  not  there  are  anomalous  data  points. 

Probability  density  calculations  of  the  white  noise  disturbance  process 
excitation  signal  n are  plotted  in  figure  49.  In  figure  49a  the  pdf  of  the  n Is 
compared  with  a theoretical  gaussian  probability  density  function  parameterized  by 
the  signal  mean  and  variance.  The  excellent  correspondence  indicates  that  n Is 
derived  from  a gaussian  process.  The  random  number  generator  from  which  n Is  di- 
rectly obtained  Is  designed  to  produce  gaussian  deviates.  The  pdf  plotted  in  fig- 
ure 49b  Is  the  same  signal  n after  It  has  been  processed  by  the  SIGANAL  TAPDC 
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Figure  49.  Signal  Probability  Density  Function  Estimate. 
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routine.  This  program  module  removes  signal  biases  and  trends  and  multiplies  the 
result  by  the  cosine  data  window  described  In  Section  III.  Multiplication  of  the 
data  by  the  window  function  alters  the  pdf  of  the  original  signal  by  adding  a con- 
centration of  new  values  near  zero;  hence  the  sharp  peak  of  the  pdf.  Otherwise 
the  pdf  retains  the  same  character  (symmetry,  mean  value)  displayed  in  figure  49a. 
Generally,  however,  pdf  calculations  should  be  performed  before  the  TAPDC  signal 
modification  required  for  transform  analysis. 

h.  Example  Summary 

The  following  conclusions  have  been  illustrated  by  the  analysis  example  and 
the  theoretical  calculations. 

• Accurate  signal  PSD's  are  obtained  from  the  magnitude  squared  of  the 
signal  discrete  Fourier  transform. 

• Averaging  must  be  employed  to  reduce  to  statistical  variation  of  PSD, 
CSD,  and  transfer  function  estimates. 

• Signal  PSD  reflects  the  signal  energy  distribution  with  frequency. 

• System  transfer  function  estimates  are  best  obtained  as  the  input-output 
signal  cross-spectral  density  magnitude  normalized  by  the  magnitude  of 
the  Input  signal  power  spectral  density. 

• The  coherence  function  indicates  the  noise  to  signal  levels  of  signals 
associated  with  a single  Input,  single  output  linear  system.  It  should 
be  calculated  to  determine  the  validity  of  transfer  function  estimates. 
The  coherence  function  also  determines  the  fraction  of  output  signal 
power  caused  by  the  Input  signal. 

• SIGANAL  Is  an  accurate  software  program  for  data  analysis  of  stochastic 
processes. 
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NOMENCLATURE 
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Ck4 
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FM 

f(y) 

F(y) 

f ,tj  ) 
f (Xj  tXy*  t^ .tg) 
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FFT 
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XX 


xy 
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(w) 

(w) 


h(t) 

H(w) 


j 


mla 


sinusoid  amplitude 
covariance £*[.(y  - uy)  (z  - uz)] 

L n 

joint  central  moment  £[(y  - u ) (z  - uz)  ] 
stochastic  process  covariance  function 
time-average  covariance  function 
Correction  factor  for  power  loss  by  data  windowing 
Discrete  Fourier  Transform 
Dirac  delta  function 

periodic  train  of  delta  functions  - perlodaT 

expectation  operator 

error  probability  confidence  bound 

cycles  per  second  frequency  variable>f=w/2ir 

Fourier  transform  of  [«] 

probability  density  function  - pdf 

probability  distribution  function 

stochastic  process  first-order  density  function 

stochastic  process  second-order  density  function 

stochastic  process  nth  order  probability  distribution 
function 

Fast  Fourier  Transform 
two-sided  power  spectral  density  - PSD 
two-sided  cross-spectral  density  - CSD 
coherence  function 

linear  dynamic  system  impulse  response 
linear  system  transfer  function,  H(w)  = F[h(t)] 
unit  Imaginary  number,  j = ^-7 
joint  moment,£[y,Kzl] 
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mxy(t,,t2) 


M 
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N 

P 

P (•) 

P(X,  f0,Af) 


*xx(f) 

*xy(f) 

Rxy(t) 


t 

T 


’d 

AT 


x(t) 


XT(t) 

X(w) 

Vk> 

(x) 


stochastic  process  correlation  function 
lower  bound  confidence  multiplier  for  PSD's 
upper  bound  confidence  multiplier  for  PSD's 
mean  value  (average)  of  the  random  variable  y 
mean  value  estimate 
number  of  data  points 
time  average  signal  power 

probability  of  the  indicated  (•)  set,  outcome  or  event 

average  power  density  of  signal  x at  frequency  f with 
bandwidth  Af  0 

phase  variable 

one-sided  power  spectral  density,  PSD 
one-sided  cross  spectral  density,  CSD 
time-average  correlation 

variance  (standard  deviation)  of  the  random  variable  y 

time  variable 

data  observation  Interval 

delay  time 

sampling  Interval 

radian  frequency  variable 

time  function  or  signal 

the  1th  sample  of  the  signal  x 

time-average  of  the  Is  ensemble  element 

truncated  time  function 

Fourier  transfer  of  x(t) 

discrete  Fourier  transform  of  x at  the  kth  frequency 
a set  or  event 
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